{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python implementation of the Kneedle algorithm\n",
    "Finding a “Kneedle” in a Haystack: Detecting Knee Points in System Behavior  \n",
    "Ville Satopa, Jeannie Albrecht, David Irwin, and Barath Raghavan  \n",
    "https://www1.icsi.berkeley.edu/~barath/papers/kneedle-simplex11.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding the knee from figure 2 from the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure2():\n",
    "    x = np.linspace(0.0, 1, 10)\n",
    "    with np.errstate(divide='ignore'):\n",
    "        return x,np.true_divide(-1, x + 0.1) + 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 0:  Raw input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = figure2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not np.array_equal(np.array(x), np.sort(x)):\n",
    "    raise ValueError('x needs to be sorted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: Fit a spline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ds = the finite set of x- and y-values that define a smooth curve, \n",
    "# one that has been fit to a smoothing spline.\n",
    "uspline = interp1d(x, y)\n",
    "Ds_y = uspline(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAa3ElEQVR4nO3deXRU553m8e+rfRegksQiQMJCIDAQY4Fj49hsTjuOE/dkbMfO4nbHHk/ScSaT5biTcU9On/Tp7tPTPU4nk22Y2EnbztJJOt326aSnByTwQgy2sI1tUCHEZlaVFrSvVfXOH1VgAQIJ1VXdulXP5xwdlW5d7v29Kunh1Vvvva+x1iIiIt6V5nYBIiISGwW5iIjHKchFRDxOQS4i4nEKchERj8tw46Q+n89WVla6cWoREc/as2dPu7W29OLtrgR5ZWUljY2NbpxaRMSzjDHHxtuuoRUREY9TkIuIeJyCXETE4xTkIiIepyAXEfE4BbmIiMcpyEVEPM6xeeTGmHSgEThprb3TqeOKiHiFtZb+kRDdg6N0D4zSMzQaeTw4Sk/0893XV7CwJN/R8zp5QdAXgSagyMFjiojEVThs6R0KvhfAY8L44lAe+7hnKEjP4CjB8OXXeDAGVi+cmZhBboypAD4M/CXwZSeOKSISi+FgiM7+ETr6RujsH6FrTPD2DI4TzkORXnTvcJArrbeTkWYozs2kODeTwtxMivOyWFCST3FuxvntRTmZ7z0e87kwO4O0NON4W53qkf898BhQeLkdjDGPAI8ALFiwwKHTikiqGBgJng/lzv4ROvpH6OwfjnzuG7st8tE3HLzssbIz0i4I2fKiHGrKC89vK8p5L5SLczMpznsvnPOy0jHG+TCORcxBboy5EwhYa/cYY9Zfbj9r7RZgC0BdXZ3WlxNJYdZa+oaD74Vv3zjhfC6wo88NjobGPVZWehqz8rOYlZ9FSUEWC0vyIo/zs5iVn31++4wxPeOczPQ4t3h6OdEjXwd81BhzB5ADFBljnrXWfsqBY4uIh/QPBznTM0Rr9xBneiIf7b3jhHP/CCPB8LjHyM5Iw1eQfT6cq0sLIo8LLgrn6LbC7IyE6yHHW8xBbq39OvB1gGiP/KsKcZHkEg5bOvpHaO0Z4kw0pFt7hjjdPXTBtt6hS4cz8rLSo73ibMqLcqidUxQN5Pd60bPys89vS8Shi0Tnym1sRSRxDI2GCPQMn+9Bn+ke5Ez3cCSgoyEd6B1iNHThiGiagdLCbGYX5bCoNJ+brimhvDiH2UWRj3OP87MVM9PN0e+wtXYHsMPJY4rI1Fhr6R4c5cy5nvOYnnSkBz3Mme5Bzg6MXvJvczPTmVOcQ3lRDmurZlFelMPsomxmF+cyOxrQvoIsMtJ1TWEi0H+VIh4WCltOdQ1ytKOfo+39HO0YiH7u58TZQYbHGYf2FWRRXpTD3OIcVi+YcUHveXY0vItyNO7sJQpykQQXDIU51TXEkY5+jnX0c6S9n2MdAxzt6Od458AFQx65meksLMljcVkhG5eWRXrSxTnne9dlhTlkZagXnWwU5CIJIBgKc+JspGd9rGMgGtaRHvbxzoELrhbMy0pnYUk+S8oL+YPls6ksyaOyJJ9KXz5lhdnqSacgBblInIyOCeujY3rVR9sjwyBjwzo/GtbL5hTxoWtnU+nLj4R1SR6lCmu5iIJcxEHBUJh3Owcu6FUf6RjgWHTMOjQmrAuyM6j05bF8XjEfXjnnfK+6siQfX0GWwlomTUEuMkWhsKUl0MfbJ7t552Q3b5/sZv+pnguuQCzMzqDSl8+KecV8dNVcFpbkU+XLY2FJPiX5CmtxhoJcZBJCYcuhtj7eOjF+aOdlpbN8bhH3r13AsrlFVPkiwyCzFNYSBwpykYucC+23T0QC+3Khfd/a+aysKGbFvGKqfAWkT8Nd7UQmQ0EuKe1qQnvFvGJWVii0JfEoyCVlXBza75zsZt+Y0M7NvDC0V8wrZlGpQlsSn4JcktJkQ/vja94bHlFoi1cpyCUpBHqHePlg+/k3I/ef7mFg5NLQXjGvmBUVxVyj0JYkoiAXT7LW4j/Ty7b9rWzzB9h7vAuIhPayuUXcW6fQltShIBfPGAmG2X2kIxLeTQFOdg0CsGr+DL5yWw0blpZRO6dIoS0pR0EuCe1s/wjbDwTY1tTKi83t9A0HyclM4+ZqH1/YWM3GpWWUFeW4XaaIqxTkknAOtfVR39TKtv0BGo91EraRBQw+smoOm5aWs67aR25Wcq25KBILBbm4LhgKs+fYWbY1tVLfFOBwez8AtXOK+PyGajbXlrNiXjFpGjIRGZeCXFzROzTKC81t1DcF2H4gQNfAKJnphvcvKuHBdZVsXFpGxcw8t8sU8QQFucTN8c6ByJBJU4DdRzoYDVlm5mWycWkZm2vL+cBiH4U5mW6XKeI5CnKZNuGwZe+JrvNDJv4zvQBcU5rPZ9ZVsXlZOasXzNQsE5EYKcjFUYMjIV5uaWfb/lbq/QHa+4ZJTzPULZzJn324lk215VT58t0uUySpKMglZoGeIbY1BahvauXllnaGg2EKszO4dUkpm2vLWb+klBl5WW6XKZK0FOQyZcc7B/hO/UF+88ZJQmHL/Fm53L92AbctK2dN5Swt8isSJwpyuWqnugb57vYWfvnacdLSDA/cuJD71iygprxAiyiIuEBBLpMW6Bni+zsO8bPd72Kx3L92AZ/fUM3sYl1ZKeImBblMqKNvmB++cIinXzlGMGy55/oKHt1YrXneIglCQS6X1TUwwpYXD/OT3x9laDTEH143jy9uWszCEs06EUkkCnK5RM/QKE++dISnXj5C30iQO1fO5YubFlNdVuB2aSIyDgW5nNc/HOQnvz/KlhcP0z04yu3LZ/Ol22pYMrvQ7dJE5AoU5MLgSIhndh3lhy8cprN/hE1Ly/jSbTVcO6/Y7dJEZBIU5ClsaDTEz199l+/vOERb7zC31JTypc2LuW7BTLdLE5GroCBPQSPBML9sPM53G1o40zPEjYtK+P4nV7OmcpbbpYnIFCjIU8hoKMxvXj/Bd+pbONk1yPULZ/LEvau4qdrndmkiEgMFeQoIhS3PvXmSb9cf5FjHAKsqivmrj63glsU+XYkpkgQU5EksHLb87p3TfGtrM4fa+lk2p4gfPVDHptoyBbhIEok5yI0x84GngXLAAlustd+O9bgyddZa/n1fK3+/rRn/mV5qygv4wSdX8wfLZ2u5NJEk5ESPPAh8xVr7ujGmENhjjNlqrd3vwLHlKlhr2X4gwBNbm3nnZA+LfPl8+773cefKuVq8QSSJxRzk1trTwOno415jTBMwD1CQx4m1lpdb2nliazNvvNvF/Fm5/N09q/jD980lI123khVJdo6OkRtjKoHrgN3jPPcI8AjAggULnDxtStt1uIMntjbz6pFO5hbn8NcfW8Hd11eQqQAXSRmOBbkxpgD4J+C/Wmt7Ln7eWrsF2AJQV1dnnTpvqhoYCfJffv4G25oClBVm8827lvPxNfPJzkh3uzQRiTNHgtwYk0kkxH9qrf2NE8eUyxsOhvjPz+xhZ0s7X/vQUh68qZKcTAW4SKpyYtaKAZ4Emqy1T8ReklzJaCjMoz97g5cOtvO3d6/knrr5bpckIi5zYiB1HfBpYKMx5s3oxx0OHFcuEgpbvvqrvWzd38o371quEBcRwJlZKy8Dmts2zay1PP7Pb/Pcm6f409uX8sCNlW6XJCIJQlMbPMBayzf/dT+/eO04X9hYzefWX+N2SSKSQBTkHvDE1mZ+vPMon1lXxZdvq3G7HBFJMAryBPeDHYf4Xw0t3LdmPv/9zlrdI0VELqEgT2BPv3KUv/m/fu5631z+8j+sUIiLyLgU5AnqV43H+cZz+7htWTl/d88q3StFRC5LQZ6AfvvWaf70n97iA4t9fPcT1+lyexG5IiVEgmnwt/LFX7zB9QtnsuXTdbrkXkQmpCBPIL9vaeezz77OsrlFPPngGnKzFOIiMjEFeYLYc6yTh59upKokn3/447UU5WS6XZKIeISCPAG8c7KbB3/8GuVFOTzz8Fpm5me5XZKIeIiC3GUHW3v59JO7KcrJ5NmHb6CsMMftkkTEYxTkLjrW0c8nf7SbjPQ0fvrwDcybket2SSLiQQpyl5zqGuQT/2c3o6EwP334Bip9+W6XJCIepSB3QaB3iE/+aDc9g6M889AN1JQXul2SiHiYo2t2ysS6BkZ44MlXOdM9xLMPr+XaecVulyQiHqcgj6PeoVH+6KlXOdzez48fXMP1C2e5XZKIJAENrcTJ4EiIh37SyL5TPXz/E6tZV+1zuyQRSRIK8jgYDoZ45JlGGo918q2Pv4/Ny8rdLklEkoiGVqbZ2MWS/8fdK/nIqrlulyQiSUY98mk0drHkP//IMu7VYskiMg0U5NPEWsuf/UtkseTHbl/Cg+uq3C5JRJKUgnwaWGv5i39t4uevHufzG67hT9ZXu12SiCQxBfk0+NbWZp7aeYQHb6rkqx9c4nY5IpLkFOQO+8GOQ3ynoYWP183nG3cu0zqbIjLtFOQOOrdY8kdWzeWvPraCNK2zKSJxoCB3yK/3nOAbz+1jc205T9yrxZJFJH4U5A747VuneezXe7m5Wosli0j8KXFidG6x5NULZrLlgevJydQ6myISXwryGJxbLLl2ThFP/fEa8rJ0oayIxJ+CfIr2HDvLw083UlmSx9Of0WLJIuIeBfkURBZLfpWywmyefegGLZYsIq5SkE/B13/zNgXZGZHFkou0WLKIuEtBfpXOdA/x9sluHrixkoqZeW6XIyKiIL9aDf4AABuXlrlciYhIhCNBboy53RhzwBjTYoz5mhPHTFQN/lbmzcilprzA7VJERAAHgtwYkw58D/gQsAy43xizLNbjJqKh0RA7WzrYVFume6iISMJwoke+Fmix1h621o4AvwDucuC4CeeVwx0MjobYoGEVEUkgTgT5POD4mK9PRLddwBjziDGm0RjT2NbW5sBp46+hKUBuZjo3LipxuxQRkfPi9mantXaLtbbOWltXWloar9M6xlpLgz/AumqfLsMXkYTiRJCfBMYuRlkR3ZZUmlv7ONk1qNkqIpJwnAjy14DFxpgqY0wWcB/wvAPHTSiadigiiSrmuzxZa4PGmEeBfwfSgaestftirizBNPhbWT63iNnFupJTRBKLI7frs9b+DvidE8dKRGf7R9hz7Cyf36BFlEUk8ejKzkl48WAbYathFRFJTArySahvClCSn8WqihlulyIicgkF+QSCoTA7DgRYv6RMiymLSEJSkE/g9Xe76BkKsqlWwyoikpgU5BOo97eSkWb4wGKf26WIiIxLQT6BhqYAa6tmUail3EQkQSnIr+B45wAHA32arSIiCU1BfgXnrubcVFvuciUiIpenIL+Cen+AKl8+Vb58t0sREbksBfll9A8H2XWoQ8MqIpLwFOSXsbOlnZFQWEEuIglPQX4Z2w8EKMjOYE3lLLdLERG5IgX5OKy11DcFuKXGR1aGvkUiktiUUuPYd6qHQO8wG5ZoWEVEEp+CfBwN/gDGwHoFuYh4gIJ8HPX+AKsqZlBamO12KSIiE1KQX6Std5i9x7s0W0VEPENBfpEdB7Q2p4h4i4L8Ig3+AOVF2SyfW+R2KSIik6IgH2MkGOalg+1sXFqGMVpEQkS8QUE+xmtHO+kbDrJxqW6SJSLeoSAfo74pQFZGGuuqS9wuRURk0hTkYzT4W7lxUQl5WRlulyIiMmkK8qjDbX0c7RjQ2pwi4jkK8qhzi0josnwR8RoFeVSDP0BNeQHzZ+W5XYqIyFVRkAM9Q6O8eqRTs1VExJMU5MBLze0Ew1bj4yLiSQpyIsMqxbmZXDd/htuliIhctZQP8lDYsuNAgPVLSslIT/lvh4h4UMon194TXXT0j+gmWSLiWSkf5Nv9AdIM3FpT6nYpIiJTkvJBXt8UoG7hLGbkZbldiojIlKR0kJ/uHmT/6R42aFhFRDwspYN8u78NQNMORcTTYgpyY8zfGmP8xpi3jDH/bIzx1Py9Bn8rFTNzWVxW4HYpIiJTFmuPfCtwrbV2JdAMfD32kuJjaDTEyy1aREJEvC+mILfW/j9rbTD65S6gIvaS4uOVwx0MjYY17VBEPM/JMfLPAP92uSeNMY8YYxqNMY1tbW0OnnZqGpoC5Gam8/5FWkRCRLxtwhUUjDHbgNnjPPW4tfa56D6PA0Hgp5c7jrV2C7AFoK6uzk6pWodYa2nwB1hX7SMnM93NUkREYjZhkFtrN1/peWPMg8CdwCZrrasBPVnNrX2c7Brk0Y3VbpciIhKzmNY0M8bcDjwG3GqtHXCmpOlX728FtIiEiCSHWMfIvwsUAluNMW8aY37oQE3Tbrs/wPK5RcwuznG7FBGRmMXUI7fWem5s4mz/CHuOneXRDZ4rXURkXCl3ZecLzW2ELWys1WpAIpIcUi7IG/wBfAVZrJxX7HYpIiKOSKkgD4bC0UUkykhL09WcIpIcUirI9xw7S89QkE26mlNEkkhKBXnDgQCZ6YabF/vcLkVExDGpFeRNAdZWzaIwJ9PtUkREHJMyQX68c4CDgT42LtVsFRFJLikT5A3+AIDudigiSSdlgrzeH2CRL58qX77bpYiIOColgrx/OMiuQx1am1NEklJKBPnOlnZGQmFNOxSRpJQSQd7gD1CYnUFd5Sy3SxERcVzSB/m5RSQ+UOMjKyPpmysiKSjpk23fqR4CvcOadigiSSvpg7y+KYAxsH5JqduliIhMi6QP8oYDAVZVzMBXkO12KSIi0yKpg7ytd5i9x7s0W0VEklpSB/n2A9GrOWsV5CKSvJI7yP0BZhflsGxOkduliIhMm6QN8pFgmBeb29iwtAxjtIiEiCSvpA3yV4900j8S0vi4iCS9pA3yBn+ArIw0bqoucbsUEZFplcRB3spN15SQl5XhdikiItMqKYP8cFsfRzsGNKwiIikhKYP83CISum2tiKSCpAzy+qYAS8oLqZiZ53YpIiLTLumCvGdolNeOduoiIBFJGUkX5C81txMMW63NKSIpI+mCvN7fyoy8TK6bP8PtUkRE4iKpgjwUtrxwoI1ba0rJSE+qpomIXFZSpd3eE1109I9oWEVEUkpSBXlDU4D0NMOtNVpEQkRSR3IFuT/A9QtmMiMvy+1SRETiJmmC/HT3IPtP92jaoYikHEeC3BjzFWOMNcb4nDjeVJy7mlOX5YtIqok5yI0x84EPAu/GXs7UbfcHqJiZS3VZgZtliIjEnRM98m8BjwHWgWNNydBoiJdb2tmkRSREJAXFFOTGmLuAk9bavZPY9xFjTKMxprGtrS2W017ilcMdDI2G2Vhb7uhxRUS8YMKbdRtjtgGzx3nqceC/ERlWmZC1dguwBaCurs7R3ntDU4DczHRuqJrl5GFFRDxhwiC31m4eb7sxZgVQBeyNDmdUAK8bY9Zaa884WuWV66PBH+DmxT5yMtPjdVoRkYQx5aEVa+3b1toya22ltbYSOAGsjmeIAzS39nGya1CzVUQkZXl+Hnm9vxXQIhIikrocW9Ay2iuPu4amANfOK6K8KMeN04uIuM7TPfKz/SO8/u5ZNi7VbBURSV2eDvIXmtsIW3S3QxFJaZ4O8np/AF9BFivnFbtdioiIazwb5MFQmBcOBNiwpIy0NF3NKSKpy7NBvufYWXqGghpWEZGU59kgb/AHyEw33LzYtRsuiogkBE8H+Q1VJRTmZLpdioiIqzwZ5Mc7BzgY6NNFQCIieDTItYiEiMh7PBnk9f4Ai3z5VPry3S5FRMR1ngvy/uEguw51aLaKiEiU54J8Z0s7I6GwFlkWEYnyXJA3+AMUZmewplKLSIiIgMeC/NwiErfUlJKZ7qnSRUSmjafScN+pHgK9wxofFxEZw1NBXt8UwBhYv6TU7VJERBKGp4J8TnEO91xfQUlBttuliIgkDMdWCIqHe9fM5941890uQ0QkoXiqRy4iIpdSkIuIeJyCXETE4xTkIiIepyAXEfE4BbmIiMcpyEVEPE5BLiLiccZaG/+TGtMGHJviP/cB7Q6W4wVqc2pQm1NDLG1eaK295B4lrgR5LIwxjdbaOrfriCe1OTWozalhOtqsoRUREY9TkIuIeJwXg3yL2wW4QG1ODWpzanC8zZ4bIxcRkQt5sUcuIiJjKMhFRDwuYYPcGHO7MeaAMabFGPO1cZ7PNsb8Y/T53caYyvhX6axJtPnLxpj9xpi3jDH1xpiFbtTppInaPGa//2iMscYYT09Vm0x7jTH3Rl/nfcaYn8W7RqdN4ud6gTFmuzHmjejP9h1u1OkkY8xTxpiAMeadyzxvjDHfiX5P3jLGrI7phNbahPsA0oFDwCIgC9gLLLtonz8Bfhh9fB/wj27XHYc2bwDyoo8/lwptju5XCLwI7ALq3K57ml/jxcAbwMzo12Vu1x2HNm8BPhd9vAw46nbdDrT7FmA18M5lnr8D+DfAAO8HdsdyvkTtka8FWqy1h621I8AvgLsu2ucu4B+ij38NbDLGmDjW6LQJ22yt3W6tHYh+uQuoiHONTpvM6wzwF8DfAEPxLG4aTKa9/wn4nrX2LIC1NhDnGp02mTZboCj6uBg4Fcf6poW19kWg8wq73AU8bSN2ATOMMXOmer5EDfJ5wPExX5+Ibht3H2ttEOgGSuJS3fSYTJvHeojI/+heNmGbo39yzrfW/jaehU2TybzGNUCNMWanMWaXMeb2uFU3PSbT5j8HPmWMOQH8DvhCfEpz1dX+vl+RpxZflghjzKeAOuBWt2uZTsaYNOAJ4EGXS4mnDCLDK+uJ/MX1ojFmhbW2y9Wqptf9wE+stf/TGHMj8Iwx5lprbdjtwrwiUXvkJ4H5Y76uiG4bdx9jTAaRP8k64lLd9JhMmzHGbAYeBz5qrR2OU23TZaI2FwLXAjuMMUeJjCU+7+E3PCfzGp8AnrfWjlprjwDNRILdqybT5oeAXwJYa18BcojcWCqZTer3fbISNchfAxYbY6qMMVlE3sx8/qJ9ngf+KPr4bqDBRt9F8KgJ22yMuQ7430RC3OtjpzBBm6213dZan7W20lpbSeR9gY9aaxvdKTdmk/m5/hcivXGMMT4iQy2H41mkwybT5neBTQDGmFoiQd4W1yrj73nggejslfcD3dba01M+mtvv7l7hXd87iPRGDgGPR7d9k8gvMkRe7F8BLcCrwCK3a45Dm7cBrcCb0Y/n3a55utt80b478PCslUm+xobIcNJ+4G3gPrdrjkOblwE7icxoeRP4oNs1O9DmnwOngVEif2U9BHwW+OyY1/l70e/J27H+XOsSfRERj0vUoRUREZkkBbmIiMcpyEVEPE5BLiLicQpyERGPU5CLiHicglxExOP+PwncicR7i3oAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, Ds_y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Normalize the spline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize(a):\n",
    "    \"\"\"return the normalized input array\"\"\"\n",
    "    return (a - min(a)) / (max(a) - min(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x and y normalized to unit square\n",
    "x_normalized = normalize(x)\n",
    "y_normalized = normalize(Ds_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Calculate the difference curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the difference curve\n",
    "y_difference = y_normalized - x_normalized\n",
    "x_difference = x_normalized.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVf7/8dcnvZAE0kBq6B0BI2CnqYAKu9hAce2sa9fd735d3fXrusX1u1/1py5WVOzYXVRcRKqoICCdhN6CQEIIAdLL+f1xbmQIgQwwyZ3yeT4e82Bm7s3M506Gd86cOfccMcaglFIq8IW5XYBSSinf0EBXSqkgoYGulFJBQgNdKaWChAa6UkoFCQ10pZQKEhroIUxEHhGRt5zrbUXkkIiE+/g5torIcB8/5hQR+atz/TwRWefLx/el473GItJcROaLyEEReUKs10SkQER+cLdyFYg00BuQE2a5IhLvcd8tIjLXxbLqZIzZboxpYoypcruWE2GM+cYY07UhHtsJ2CdEJN+5fHgqj1fHazwR2AskGmN+C5wLXAi0NsYMOLXqVSjSQG944cA9p/ogTrjo76txXQRMAE4HWgIv+vjx2wFrzeGz+9oBW40xRSf6QCIS4dPKXBZsx9NYNCAa3j+B34lI07o2isjZIrJYRAqdf8/22DZXRP4mIt8CxUAH576/ish3zsf3z0QkRUTeFpEDzmNkeDzG0yKyw9m2VETOO0YdGSJiRCRCRM5yHrvmUioiW539wkTkARHZ5LRa3xeRZI/HuU5EtjnbHjreCyMio0RkrdPlsFNEfufcP1hEckTkQRHZ63zSufYYjzFYRHI8bm8Vkd+JyErnNX1PRGI8tl8qIstFZL/zGvY5TokVQAmw2xhTZoyZebzjcR6/vYjMc45pJpDqsc3zNZ4CXA/83nmNfw1MBmpe+z/XV69zrP8tIiuBIudxW4rIRyKSJyJbRORuj/0fcX5fbzj1rRGRTI/tbUTkY+dn80XkXx7bbhKRLLHdQTNEpN1xXoNznVr3O++9G5z754rILR773SAiCzxuGxG5Q0Q2ABtE5HkR+b9aj/1vEbnfuX7MYw1Zxhi9NNAF2AoMBz4G/urcdwsw17meDBQA1wERwHjndoqzfS6wHejpbI907tsIdASSgLXAeud5IoA3gNc8apgApDjbfgvsBmKcbY8AbznXMwADRNQ6hkhgHvCYc/seYCHQGojGtlrfdbb1AA4B5zvbngQqgeHHeH12Aec515sB/Z3rg52fe9J5nAuAIqCrs32Kx+s5GMip9Zr/gG1RJwNZwG3Otn5ALjAQ+8npemf/6GPU1xI44DxfmJe/8+896j4fOHis19jzOJzbNwALPG4ft17n+nKgDRCLbaAtBR4GooAOwGbgYo/fdykwynm8x4CFzrZwYAXwFBAPxADnOtvGYN9z3bHvoz8C3x3j+Ns5xzwe+95JAfp6vJ9vOc7xGmCm83uLdV6/HYB4vEdKnN/LcY81VC/aQm8cDwN3iUharfsvATYYY940xlQaY94FsoHLPPaZYoxZ42yvcO57zRizyRhTCHwJbDLGfG2MqQQ+wAYBAMaYt4wx+c7PP4ENmhPpc34G+x+0prV9G/CQMSbHGFOGDYkrxH5EvgL43Bgz39n2J6D6OI9dAfQQkURjTIEx5sda2/9kbMt4HvAFcJW3NRtjfjLG7AM+A/o6908EXjTGLDLGVBljXgfKgEG1H0BEIoEZwO3YIJksTpeXiCwQkcvq+Jm2wJkedc93nv9keVPvM8aYHcaYEue504wxjxpjyo0xm4GXgXEe+y8wxkw3th//TWx3EsAAbFD+lzGmyBhTaoypaT3fhv2DnuW8x/4O9D1GK/0a4GtjzLvGmArnvbf8BI75MWPMPud4vsGGfM2nyiuA740xP3l5rCFHA70RGGNWA58DD9Ta1BLYVuu+bUArj9s76njIPR7XS+q43aTmhtP9kOV0P+zHtupT8YLTDTAYuMYYUxPM7YBPnI/T+7Et4CqguXM8P9drbF9w/nGe4nJsa3Gb001xlse2AnNkX/I25/G9sdvjejGHX492wG9ranfqb3OMxx0KRBlj3gKuBtpjQz0R6AYsqONnWh6j7pPlTb07au3fstb+D2J/NzVqvzYxzh/jNsA2J7DrquNpj8fcBwhHvk9rtAE2ndhhHsHz/WOAqdjWPtg/Fm971FTfsYYc/eKh8fwP8CPwhMd9P2HfmJ7aAv/xuH3S02GK7S//PTAMWGOMqRaRAux/Rm9+9i/Yj90HPDbtAG4yxnxbx8/swn4sr7kdh/3IXSdjzGJgjNMavhN4HxsIAM1EJN4jHNsCq+urux47gL8ZY/7mxb41XVwYY0pFZDQwB1gMTDXGFNTxM7uOUffJ/g69qdfzsXcAW4wxnU/yudqKSEQdoV5Tx9t1/Fxdj3OsETpFQJzH7RZ17FP7tXoX+EpE/oHtevqlx/Oc7LEGLW2hNxJjzEbgPcDzi5vpQBcRucb5QutqbD/05z562gRsX3QeECEiDwOJ9f2QiLTBhuuvjDHra21+AfhbzcdtEUkTkTHOtg+BS50vxaKARznGe0xEokTkWhFJcrqSDnB098yfnf3OAy7FdiedipeB20RkoFjxInKJiCTUse8CbOv1URGp6Z+eA3TBtmyPYozZBizxqPtcjuw+a8h6wX53cFDsF6WxIhIuIr1E5EwvnusH7B+kfzjPEyMi5zjbXgD+ICI9AUQkSUSuPMbjvA0MF5GrnPd0iojUdHktB8aKSJyIdAJurq8oY8wy7NDOycAMY8x+Hxxr0NJAb1yPYr9wAsAYk48Nqt9iuyZ+D1xqjNnro+ebgW3tr8d+9C+l7i6c2oZhP7p+KIdHuqxxtj0NTMO2mg5ivyAd6BzPGuAO4B1sOBQAObUf3MN1wFYROYDtp/UcybLb+fmfsCFxmzEm24vaj8kYswS4FfiX89gbsV/M1bVvIXbY4iCnhk3YTxsDgBtF5NZjPM012NdjH/ZT2RuNUa+zfxX2/dQX2MLhIEzy4rmqsH98OmG/iM/BdjVhjPkEeByY6vyuVgMjj/E427HdaL/FvgbLOdxP/xRQju0ifJ3D3Sf1eQf7pf87vjjWYFbz7bFSfkNEBmNHhrR2uxalAom20JVSKkhooCulVJDQLhellAoS2kJXSqkg4do49NTUVJORkeHW0yulVEBaunTpXmNM7bPOARcDPSMjgyVLlrj19EopFZBE5JhnH2uXi1JKBQkNdKWUChIa6EopFSQ00JVSKkhooCulVJCoN9BF5FWxCx3XOXWpMwvcMyKyUeyyX/19X6ZSSqn6eNNCnwKMOM72kUBn5zIReP7Uy1JKKXWi6h2HboyZLx6LDtdhDPCGs7rIQhFpKiKnGWN2+ahGpZQKGMYYSiqqKCyp4EBJJYUlFUddhnVL5/Q2da4bf0p8cWJRK46cYzvHue+oQBeRidhWPG3btvXBUyullO8ZYzhYVklhcQUHSiuccD4cyHUF9YHSw/tUVB1/jqz0hGi/DXSvGWNeAl4CyMzM1FnBlFINqryymoLicvIPlbOvqJz9JeVHhfKBWqFcc1/1cRIqPExIjIkgKTaSxNhIkmIjadUsliTnemJM5M/Xf74v1u6fEBNJeFi9q0CeFF8E+k4OrwMJ0Nq5TymlfKq0oop9RTac84vK2VdU9nNY19yXf6js5+sHS+ta89qKCg8j0SNoU5pE0SEt/qhATqwVyEmxkTSJjkCkYUL5VPgi0KcBd4rIVOzSW4Xaf66U8kZxeeVRgbyvqMwJ5iPv23eonKLyqjofJyJMaBYfRUp8FMnxUfRqleRcjyalyeH7m8ZF/RzKMZFhfhnKp6LeQBeRd4HBQKqI5GDXSaxZDf0F7ELHo7DrHRYDNzZUsUop/1daUcWeA6XsLixlt/Nv3sEyj3CuCeoySitqrwtuRYWHkeyEcEqTKDJS4uz1OkI6JT6axFj/bDE3Nm9GuYyvZ7vBLgyslApixhj2F1f8HNI1/+45UHrEffuLK4762ZjIMFLio38O6c7pTez1JnWHtL92afg716bPVUr5j/LKanIPOuFcWMbuA/b6rsJS9jhBvedAKWWVR7aoRSAlPpoWSdG0bhbLGe2a0SIxhhZJziUxhuZJMSRoQDcKDXSlgljN8Ls9Hi3q3R4BvdsJ8PyiMmqvRhkVEWbDOTGGvm2a0iIphubO7RZJ0bRIiiU9IZrIcJ1BxF9ooCsV4KqrDbsPlLI1v4ht+cVs3VvE1vwitu4tJqeguM4vEpvFRdpwToqhV8skmifGcFqSbU3XhHjTuEhtVQcYDXSlAkB1tWHXgVK27S1iS63g3pZffERXSFR4GG1T4shIiePsTik2qH9uWdvrMZHhLh6Naiga6Er5iapqw67CErbuLXaCuogte4vZll/Etn3FlHuGdkQYGSlxtEuJ54IuaWSkxpOREk+7lDhOS4ptsBNXlH/TQFeqEVVVG37aX2K7RJxW9jbn+vb8YsqrDod2dEQYGSnxtE+NZ2i3dNqlxJOREkdGajwtEmMI09BWtWigK+Vj1dWGnftL2LL3yFb2lvwiduwrPmKej5hIG9qd0powrHs6GSm2pZ2RGkfzBA1tdWI00JU6BdXVhm37ilm1s5BVOftZtbOQNTsPcLDs8CnnsZHhtEuJo2vzBC7q0YL2qXFOazue5onR+sWj8hkNdKW85Bneq3cWsjJn/xHhHRURRvfTEhnTryU9WybRIdV2l6QlaGirxqGBrlQdaof3qhz7b13h3btVEr1aJdGleYKOyVau0kBXIa/O8P6p8OeZ+qIiwujeIkHDW/k9DXQVUowxbMsvZuWxwjs8jO6nJTD69Jb0aa3hrQKLBroKWjXhvWpnofOl5bHD27PlHRWh4a0Ckwa6Chr7i8v5dmM+K3P2s7KO8O6m4a2CnAa6Cmhb9hYxK2sPM9fuYcm2Aqqqzc/hfdnpLemj4a1CiAa6CihV1YYftxfwddYevl67h015RQB0bZ7Ar8/vwLDu6fRu1VTDW4UkDXTl9w6VVfLN+jxmZu1h7ro89hWVExEmDOqQwoRB7RjevTltkuPcLlMp12mgK7+0c38Js7L28HVWLgs35VNeVU1SbCRDuqYxvEdzzu+SRmJMpNtlKuVXNNCVX6iuNqzaWWj7w7Nyydp1AID2qfFcf3Y7hnVvTma7ZkTo8EGljkkDXbmmtKKKbzfu5eusXGZl7SH3YBlhApntkvnDyG4M79GcjmlN3C5TqYChga4aVe7BUuZk5zJzbS4LNuZRWlFNfFQ4F3RNY3j35gzumk5yfJTbZSoVkDTQVYMyxrBuz0G+Xmv7w5fv2A9Aq6axXJ3ZhmHdmzOwQzLREbqCjlKnSgNd+Vx5ZTWLtuQzKyuXmWv3sHN/CQCnt2nKby/swvAezenWIkFnIFTKxzTQlU/sLy5nzrpcvs7KZf66PA6WVRITGca5nVK5a2gnhnZLJz0xxu0ylQpqGujqlOw9VMYLczfx5sJtlFVWk5YQzaWnn8awbs05p1MqsVHalaJUY9FAVyeloKicF+dv5vXvtlJWWcXY/q2ZMKgdfVol6bJpSrlEA12dkMKSCl75ZjOvfruVovJKRp/eknuGdaaDDi9UynUa6Morh8oqeW3BFl7+ZjMHSisZ1bsF9w7vQpfmCW6XppRyaKCr4your+SN77fx4rxNFBRXMLx7c+67sDM9Wya5XZpSqhYNdFWn0ooq3l60nefnbmLvoTIu6JLG/Rd24fQ2Td0uTSl1DBro6ghllVW8v3gH/5qzkT0Hyji7YwovTOhPZkay26UpperhVaCLyAjgaSAcmGyM+Uet7W2B14Gmzj4PGGOm+7hW1YAqqqr5aGkOz87eyM79JZyZ0Yz/d3U/zuqY4nZpSikv1RvoIhIOTAIuBHKAxSIyzRiz1mO3PwLvG2OeF5EewHQgowHqVT5WVW34dNlOnpm9gW35xZzepimPje3NeZ1T9UxOpQKMNy30AcBGY8xmABGZCowBPAPdAInO9STgJ18WqXyvutrw+apd/L+v17M5r4ieLRN55fpMhnZL1yBXKkB5E+itgB0et3OAgbX2eQT4SkTuAuKB4XU9kIhMBCYCtG3b9kRrVT5gjGHGmt08NXMD6/YcpGvzBF6Y0J+Le7bQIFcqwPnqS9HxwBRjzBMichbwpoj0MsZUe+5kjHkJeAkgMzPT+Oi5lReMMczOzuXJmetZ89MBOqTF88z4flza+zQ9s1OpIOFNoO8E2njcbu3c5+lmYASAMeZ7EYkBUoFcXxSpTp4xhm827OXJmetZvmM/bZPjeOLK0xnTt6Wu/qNUkPEm0BcDnUWkPTbIxwHX1NpnOzAMmCIi3YEYIM+XhaoT9/2mfJ6cuY7FWwto1TSWf4ztzeVntCZSg1ypoFRvoBtjKkXkTmAGdkjiq8aYNSLyKLDEGDMN+C3wsojch/2C9AZjjHapuGTptn088dV6vtuUT/PEaP7yi15cndmGqAgNcqWCmVd96M6Y8um17nvY4/pa4BzflqZO1Iod+3ly5nrmrc8jtUkUf7q0B9cObEtMpE5hq1Qo0DNFg8Canwp5auYGvs7aQ7O4SP4wshvXndWOuCj99SoVSvR/fACrqjY89Mkqpi7eQWJMBL+7qAs3nNOeJtH6a1UqFOn//ABVXW34749W8uHSHG49rz13Du1MUmyk22UppVykgR6AjDH8+bM1fLg0h/uGd+Ge4Z3dLkkp5Qd02EOAMcbw+H/W8fr32/j1+R24e1gnt0tSSvkJDfQAM2nORl6Yt4kJg9rywMhuerq+UupnGugB5JUFW/i/r9Yztl8rHh3dS8NcKXUEDfQAMfWH7fzl87WM7NWC/72ij86/opQ6igZ6APj38p384ZNVDO6axtPj+ukcLEqpOmky+Lmv1uzm/vdXMLB9Mi9MOENP31dKHZOmgx+bvz6PO99ZRu9WSUy+/kw9hV8pdVwa6H7qhy37mPjmEjqmN+H1Gwfo2Z9KqXppoPuhlTn7uWnKYlo1jeXNmweQFKdngCql6qeB7meydx/gV6/+QLP4SN6+ZRCpTaLdLkkpFSA00P3I5rxDTJj8AzER4bxzyyBaJMW4XZJSKoBooPuJHfuKuXbyIowxvHXLQNokx7ldklIqwGig+4HcA6VMeGURRWWVvHnzQDqlN3G7JKVUANJAd9m+onKunbyIvQfLmHLTAHq0THS7JKVUgNKxcC4qLKngulcWsX1fMVNuHED/ts3cLkkpFcC0he6SorJKbpqymPV7DvLCdWdwVscUt0tSSgU4DXQXlFZUcesbS1i2vYBnxvVjSNd0t0tSSgUB7XJpZBVV1dzx9o98tymfJ686nZG9T3O7JKVUkNAWeiOqqjbc+95yZmXn8tdf9GJs/9Zul6SUCiIa6I2kZlHnL1bu4qFR3ZkwqJ3bJSmlgowGeiPwXNT5nmGdufX8Dm6XpJQKQhrojeB/Z9hFnW89rz33Du/sdjlKqSClgd7AJs3ZyPNzN3HtwLY8OKq7rgOqlGowGugN6NUFW/jnjHX8sl8r/jJGF3VWSjUsDfQGMvWH7Tz6+Vou7tmcf+qizkqpRqCB3gBqFnW+oEsaz4zXRZ2VUo3Dq6QRkREisk5ENorIA8fY5yoRWSsia0TkHd+WGThqFnUekGEXdY6O0HVAlVKNo94zRUUkHJgEXAjkAItFZJoxZq3HPp2BPwDnGGMKRCQkz2X/ZoNd1LlXqyReueFMYqM0zJVSjcebFvoAYKMxZrMxphyYCoyptc+twCRjTAGAMSbXt2X6vx+27OPWN5bQIS2e1288Uxd1Vko1Om8CvRWww+N2jnOfpy5AFxH5VkQWisiIuh5IRCaKyBIRWZKXl3dyFfuhmkWdWzaN5c2bB9I0LsrtkpRSIchX39ZFAJ2BwcB44GURaVp7J2PMS8aYTGNMZlpamo+e2l3rdh/kV6/+QNO4SN6+ZSBpCbqos1LKHd4E+k6gjcft1s59nnKAacaYCmPMFmA9NuCD2pa9RVw7eRHREWG8c8sgTkuKdbskpVQI8ybQFwOdRaS9iEQB44Bptfb5FNs6R0RSsV0wm31Yp9/JKSjm2pcXUm0Mb98ykLYpuqizUspd9Qa6MaYSuBOYAWQB7xtj1ojIoyIy2tltBpAvImuBOcB/GWPyG6pof/C3L7I4WFrJGzcNoFN6gtvlKKWUdwtcGGOmA9Nr3fewx3UD3O9cgl5ZZRXz1ufxy36t6NUqye1ylFIK0DNFT8qizfsoLq9iaLeQHG6vlPJTGugnYXZ2LtERYZzdMdXtUpRS6mca6CfIGMOs7D2c0ylVzwRVSvkVDfQTtCnvEDv2lTBEu1uUUn5GA/0Ezcqysxpo/7lSyt9ooJ+g2dm5dGuRQKumehKRUsq/6AxSJ6CwuIIl2wr4tT8s8lxRClnTYNmbUHoAkttDs/aQ3OHw9YTTIEz/ZisVKjTQT8D8DXlUVRuGdXexuyVvPfz4Oix/B0r2QbMMG+I/LYe108BUHd43IsZurx30ye2haVsIj3TrKJRSDUAD/QTMzs6lWVwkfds0a9wnriyzYb10CmxbAGER0O1SOOMGaH/B4VZ4VQUU7oB9W6Bgi/235vrmuVBZcvgxJRySWttwT+5wOOiTO9g/AlHxjXuMSqlTpoHupapqw9x1uQzumk54Y60PuneDDXHP1viw/4F+E6BJHZ8SwiOdlngdXULGwKE9sG9zrcDfDGs+gZKCI/dv0ryOlr1zPbYZ6ILXSvkdDXQvLd9RQEFxRcOPbqksg6zPbJBv/cZpjV/itMYHn3yfuAgktLCXdmcfvb2k4Migr/l381xYUWtFwegkpzXf/siWfcv+EKWTlCnlFg10L83KyiU8TDi/SwPN4753I/w4xbbGi/OhaTvbGu97LSQ0b5jn9BTbDFo1g1b9j95WUQIF22xr3jPwd62wf3yqK+1+kfHQbRT0ugI6DoUIXehDqcakge6l2dm5ZLZrRlKsD79IrKs13nUUZN54aq1xX4uMhfRu9lJbVSUcyLFf1mZ/Dmv/Das+gJim0GO0DfeMcyFMz6pVqqFpoHth5/4Ssncf5MFRdQTaycjfBEtfq9Uafxj6Tmic1rgvhUc4I2kyoMtFMOr/YPMcWPUhrPoIfnzD9sf3/KUN99aZ2v+uVAPRQPfC7Oyas0NPIWwryyHbaY1vmW9HmXQbBWfcCB2G+E9r/FRFREGXi+2lvBg2zLDhvuQ1WPSC/ePV63J7ad5Tw10pH9JA98LsrD20TY6jY9pJDOXL32THjS97G4r32vHfQ/9kR6oktPB9sf4kKs62zHv+EkoLIfsLG+7fPg0LnoS0brbV3msspHR0u1qlAp4Gej1Kyqv4blM+4we0RbxtTVaW2/7kpVNgyzzbGu860vaNdxgaPK3xExGTBH2vsZdDebD2U1j9Ecz5q7207Hc43BNbul2tUgFJA70e323aS1lltXdnh9ZujSe1haF/hH7XBX9r/EQ0SYMBt9pLYQ6s/hhWfwhfPQRf/RHanWODvccvID7F7WqVChga6PWYlZ1LXFQ4A9on171DZTms+8K2xjfPPdwaP+NG6DhER3fUJ6k1nHO3vezdaFvtqz+EL+6HL39vv1/odbkdix+T6Ha1Svk1DfTjMMYwJzuX8zqnEh1RK5iL9sL3/4Jlb0FRHiS1gSF/tH3jiae5U3CgS+0Eg/8bLvg97Flt+9tXfwyf3mbnpel8EfS+wv4bqbNdKlWbBvpxZO06yK7CUu4d3vnojZ/eDhu/dlrjN9gTabQ17hsi0KK3vQx/BHb8YFvuaz6xM0xGJdgWe+8roMNgnWRMKYcG+nHMWWeHKw7pWqv//OAe2DgTzr3Pjh9XDUcE2g60l4v/bk/AWv0hrP0MVk6F2GToMcaGe9uzQ/MLZ6UcGujHMStrD31aJ5GeGHPkhlUfgKmGPuPcKSxUhUfY7yU6DoFLnoSNs2y4r3zPnqiV0NJ+mZp5kw6DVCFJA/0Y8g+VsWzHfu4eWkd3y4qp0OoMSOvS+IUpKyLanpjVbRSUF8G6L223zKIX4ftJdgqFs++EtmfpyUsqZGigH8O89XkYw9HDFXevhj2rYOQ/3SlMHS0q3na59L4CDuXCDy/D4sl29FHL/jbYu4+xLXylgph2OB7DrOxc0hKi6dUy6cgNK6faSbR6Xe5OYer4mqTD0IfgvjVwyRP2DNUPb4Jn+tmWe9lBtytUqsFooNehoqqa+evyGNI1jTDPxSyqq2DlB3bYnJ7w4t+i4uDMW+DOJTDuHTvefcaD8GQPe/JSYY7bFSrlcxrodViytYCDZZVHT8a1eS4c2g2n65ehASMszA5xvOlLuHU2dBoO3z8HT58OH91i12JVKkhop2IdZmfvISo8jHM7px65YcVUOydJlxHuFKZOTasz4MrXYP92WPiCndp31QeQcR6cdaf95KXDHlUA03dvHWZl5zKwQzJNoj3+3pUdtBNu9RxrR1iowNW0LYz4O9y/Bi78i12J6d2r4bmBdprfipL6H0MpP6SBXsvWvUVszis6eu3QrM+gohhOH+9OYcr3YpLsHDL3rICxk+10Ap/fC0/1gjmP2VkhlQogXgW6iIwQkXUislFEHjjOfpeLiBGRTN+V2LgOL2ZRK9BXTLULIrcZ4EJVqkGFR0KfK2HiPLj+c7uq0rx/wFM9Ydrddnk9pQJAvYEuIuHAJGAk0AMYLyI96tgvAbgHWOTrIhvT7OxcOqbF0y7FYzGLwhy7ytDp4/QklWAmAu3Pg2vegzsWQ9/x9izUSWfC21fZ94Axblep1DF500IfAGw0xmw2xpQDU4Exdez3F+BxoNSH9TWqQ2WVLNqSz7DutUa3rHwfMNDnalfqUi5I6wKXPW3Hsw/+A+xcCq9fBi9dYN8PVRVuV6jUUbwJ9FbADo/bOc59PxOR/kAbY8wXx3sgEZkoIktEZElenv/1Ty7YkEdFlTmyu8UY20prMwiS27tXnHJHfCoMfgDuW20DvqIEPr7VDnv89ml74pJSfuKUvxQVkTDgSeC39e1rjHnJGJNpjMlMS0s71af2udnZuSTERHBGu2aH79y1HPKydex5qIuMtdMk374IrnkfkjvAzIfhyZ7wnwftUEilXG4YNqgAABRISURBVOZNoO8E2njcbu3cVyMB6AXMFZGtwCBgWqB9MVpdbZidnccFXdKIDPd4WVa8B+HR0PMX7hWn/EdYGHS5GG743H6J2nUk/PAiPN0XPrgBcpa6XaEKYd4E+mKgs4i0F5EoYBwwrWajMabQGJNqjMkwxmQAC4HRxpglDVJxA1m1s5C9h8qOnIyrqsKeeNJ1BMQ2O/YPq9DUsi9c/rId9njWHXY638lD4dURkP0FVFe7XaEKMfUGujGmErgTmAFkAe8bY9aIyKMiMrqhC2wss7NzEYELungE+sZZdrFnHXuujiepNVz0F7h/LVz8GBTuhKnXwEvnQ/Z0HRmjGo1Xp/4bY6YD02vdV+dSPcaYwadeVuObnZ1L/7bNSI6POnznyqkQl2Ln/1CqPtEJcNbtMGCiXXhj3uMwdTyc1heGPGinFtBhr6oB6ZmiQO6BUlbtLDxydEvJftu66nWFrlmpTkx4hP0S/Y7FMOY5KCmAd66CycPsOrTaYlcNRAOdw2uHHhHoaz+FqjId3aJOXngE9LsW7loKlz1jF99463J49WI7c6cGu/IxDXRgVlYuLZNi6NYi4fCdK96D1C7Qsp97hangEB4JZ1wPd/1o10ItzIE3xsCUS2DrArerU0Ek5AO9rLKKBRv3MqRbOlLTv1mwFbZ/p6f6K9+KiIIzb7bBPvKfkL/Jhvrrl8G2792uTgWBkA/0RZv3UVxedeRwxZXv2397X+VOUSq4RcbAwIlwz3IY8Q/IzYbXRsAbv4AdP7hdnQpgIR/os7NziYkM4+yOzmIWxsCKd+2iB03bHP+HlToVkbEw6Dd2HPtFf4Xdq+CVC+GtK+zcMUqdoJAOdGMMs7L3cHbHVGIiw+2dOYvtggc69lw1lqg4OPsuG+zDH4GdS+DlofDOONi1wu3qVAAJ6UDflHeIHftKjhzdsmIqRMRCj6A5Z0oFiugmcO59cM9KGPpH+z3Oi+fD1Gth92q3q1MBIKQDfVZWreGKlWWw+iPofqk9SUQpN8Qkwvn/BfeuslP3bpkPL5wD718PuVluV6f8WEgH+uzsXLq1SKBl01h7x4avoHS/jj1X/iEmyU7de+9KOP/3diqK586CD2/SVZRUnUI20AuLK1iyreDI0S0rpkKT5tB+sGt1KXWU2GYw9CEb7OfeB+v+Yxe0/niiHfqolCNkA33ehjyqqg1DuzmrExXvg/UzoPeV9gw/pfxNXDIM/x8b7GfdCWunwb/OhE9vh31b3K5O+YGQDfQ52bkkx0fRt01Te8fqj6C6QrtblP+LT7WzO967EgbeZt+7z54B0+6Cgm1uV6dcFJKBXlVtmLMul8Fd0ggPc84EXTEVmveCFr3dLU4pbzVJhxF/t8MdB9xq38PPngGf3WunF1AhJyQDfdn2AvYXVzC0pv9870Y79ldb5yoQJbSAkY/D3cvtnDHL3oJn+sEXv4MDP7ldnWpEIRnos7NzCQ8TzuvsrGu6cipImO0/VypQJbWCS56Au5dB32tg6Ws22Gf+j50OWgW9kA30MzOakRQbaZcJW/EedBhiWzpKBbqmbeCyp+20vT3GwLdPwzN94ftJ9lwLFbRCLtB37i8he/dBhtWMbtn+PRRu1+4WFXyaZcDYl+DX8+000DMehH9l2snndL3ToBRygT47254dOqTm7NAV70JUE+h2iYtVKdWATusD130C130KMU3h41vhpQtg02y3K1M+FnqBnrWHdilxdEyLh4oSWPtv+7E0Kt7t0pRqWB2HwMR5MHayPSP6zV/aKXt1ArCgEVKBXlJexXeb8hlas5jFuulQdgD6XO12aUo1jrAw6HMl3LkELn4Mdi23E4B9dKuOYQ8CIRXo323aS1ll9eHJuFZMhcTWdu5zpUJJRDScdbsd6njufZA1zfav/+dBe9a0CkghFeizsnOJjwpnQPtku2DvxlnQ5yrbalEqFMU2tXOw3/Wj/b+w6Hl4ui9886TtklQBJWSSzBjDnOxczu2cSnREOKz6EEyVjm5RCuwY9jGT4DffQbuzYNaf4Zn+8OObUF3ldnXKSyET6Fm7DrKrsPTwcMUV79qhXGld3S1MKX+S3h2ueQ9umA6Jp8G0O+H5c+zEdca4XZ2qR8gE+uzsPQAM7pYGe9bC7pXQR1vnStUp4xy4ZRZc+TpUlcE7V8GUSyFH1zr1ZyEU6Ln0aZ1EekKMPdU/LAJ6Xe52WUr5LxHo+Qu44wcY9X+Qlw2Th9qVk3Qedr8UEoGef6iMZTv229Et1VWw8gPodCE0SXO7NKX8X3iknc3xnuVwwQOwYSZMGmAn/zqU53Z1ykNIBPrcdXkYg+0/3zIfDv4Ep+vYc6VOSHQCDPmDnfyr//Ww5FU7R8zcx6HskNvVKUIk0GevyyUtIZqeLRPt2PPoJOgy0u2ylApMCc3h0ifhjkX27NO5f7ezOi5+Baoq3K4upHkV6CIyQkTWichGEXmgju33i8haEVkpIrNEpJ3vSz05FVXVzF+Xx9Cu6YRVFEHWZ7ZfMDLG7dKUCmypneHqt+DmmZDSEb64H54bZJfG0xExrqg30EUkHJgEjAR6AONFpEet3ZYBmcaYPsCHwP/6utCTtXjrPg6WVdrFLLI/h4oiOH2822UpFTzaDIAbv4Rx74KEw/vXwSsXwbbv3a4s5HjTQh8AbDTGbDbGlANTgTGeOxhj5hhjip2bC4HWvi3z5M3JziUqPIxzO6Xa7pam7aDtILfLUiq4iEC3UfbEpMuegcId8NoIePcayFvndnUhw5tAbwXs8Lid49x3LDcDX9a1QUQmisgSEVmSl9c4347Pys5lYIdk4styYfNce2aoSKM8t1IhJzzCLoN3148w9E92EMJzg+Cze+x0G6pB+fRLURGZAGQC/6xruzHmJWNMpjEmMy2t4YcMbt1bxOa8IoZ1S7eT+mN0ZkWlGkNUHJz/O2cB64nOOqf9YcFTUFHqdnVBy5tA3wm08bjd2rnvCCIyHHgIGG2M8Yt1rmoWsxjaNd12t7QZaL+8UUo1jvgUu4D17Qsh41z4+hGYdCas+US/OG0A3gT6YqCziLQXkShgHDDNcwcR6Qe8iA1zv/lcNTs7l07pTWhbsRHysrR1rpRbUjvDNVPhV/+GqAT44AZ4bSTs/NHtyoJKvYFujKkE7gRmAFnA+8aYNSLyqIiMdnb7J9AE+EBElovItGM8XKM5VFbJoi35trtlxVQIj4Kev3S7LKVCW4fBcNs3dhHr/I3w8hD45DY48JPblQWFCG92MsZMB6bXuu9hj+vDfVzXKVuwIY+KKsOQLsnwyQfQ5WKIS3a7LKVUWDiccQP0HAvfPAELn7NLQZ5zD5x9t+1/VyclaM8UnZWVS2JMBJmVy6AoT8eeK+VvYhLhwj/DnYuh80Uw9zF49gxY8R5UV7tdXUAKykCvrjbMWZfHBV3TiVj9PsQm28m4lFL+p1kGXPW6PTmpSTp8MhEmD4PtC92uLOAEZaCv2lnI3kNlXNQxBrK/sNPkRkS5XZZS6njanQ23zoFfvAAHd8GrF9svT3Xxaq8FZaDPys4lTGBI1fdQWardLUoFirAw6Dse7lpqp+pd9x/415nw9Z+h7KDb1fm9oAz0Odm59GvbjCbZH0FKJ2jV3+2SlFInIireTtV711I7md6CJ+2JSUtf1zVOjyPoAn3PgVJW7SxkTEYlbFugp/orFciSWsHYl+CW2ZDcHj67G168ADbPc7syvxR0gT7HOTt0ZPV8e4eeTKRU4Gt9Btw0A654DUoL4Y3RduIvXQrvCEEX6LOzc2mZGE3q5k+g3bnQtK3bJSmlfEEEeo21wxyHPQxb5sGkgfCfB6GkwO3q/EJQBXppRRULNu7lV+32IvkbbXeLUiq4RMbAeb+1Mzr2HW9PTHqmPyx6KeRXTAqqQF+0ZR/F5VVcyjcQEQM9xtT/Q0qpwJTQHEY/a6cSaN4TvvwveP4cu4h1iAqqQJ+TnUtCZDWtcr6AbpfYM9GUUsGtRW+4/jMY9w5UV8DbV8CbYyE3y+3KGl3QBLoxhlnZe5h42iakpEDHnisVSkRsI+72RXDx32HnEtta//x+KNrrdnWNJmgCfVPeIXbsK2GMzIf4dOgwxO2SlFKNLSIKzroD7loGZ94MS6fY/vXvnoXKcrera3BBE+izsnJJ4hBt8uZD7yvtUlhKqdAUnwKj/mnXOG0zAL76I0waYKcCCeKFNYIn0LNzuSV5OVJdoaNblFJWejeY8CFM+AgiomHqNfDmL4N24eqgCPTC4gqWbitgbNg3kN7DfkmilFI1Og2H2xbAiMfhpx/h+bPhP3+Akv1uV+ZTQRHo8zbk0drsotWhVXqqv1KqbuGRMOg2O3693wRY+Lydfz2I5ocJikCfnbWHa2O+xyC2/1wppY4lPtUugffreXat08/utkvhBcH86wEf6FXVhrnrcrk8YgHSYTAktnS5IqVUQDjtdLuoxuWvwKE8O//6R7cG9PqmAR/oy7YX0Kl0NSkVu3TsuVLqxIhA7yvgriVw3u/s2qbPZtq1TitK3a7uhAV8oM/KzuWKiAWYyHjofqnb5SilAlFUPAz7E9yxCDoOgVmPwnODIHt6QA1zDPhAX5CVw2URi5Aeo+0vRSmlTlZyexj3Nlz3CYRHwdTx8NbYgBnmGNCBnlNQTJu8ecSbIp33XCnlOx2Hwm++hRH/gJylzjDHB+1c7H4soAN9TnYuY8O/oTK+BbQ/3+1ylFLBJDwSBv0G7v4R+l5rp+l99gz48Q2orna7ujoFdKAvXrOOweErCe87DsLC3S5HKRWM4lNh9DMwcQ4kd4Bpd8HkobDjB7crO0rABnpJeRXp274ggipET/VXSjW0lv3sMnhjJ8PB3fDKhfDxr+HALrcr+1nABvp3m/YyWuZzKLknpHd3uxylVCgQgT5Xwp1L7KpJaz623TALnoLKMrerC9xAX7V8EX3CthCTOcHtUpRSoSa6iV3X9I5F0GEwfP2IHea47ktXhzkGZKAbY0je+AlVhBHRR0/1V0q5JLkDjH8HJnwMYRHw7ji7YtLeDa6UE5CBnvVTIcMr57In/VxokuZ2OUqpUNdpmJ17/eLH7Jelzw2CGQ81+jDHgAz0dQun01L2EX/mtW6XopRSVngknHW7nc2x7zXw/STbv77srUYb5uhVoIvICBFZJyIbReSBOrZHi8h7zvZFIpLh60I9JW34iCKJI6nvmIZ8GqWUOnFN0mD0s3aYY7P28O87YPIw2LG4wZ+63kAXkXBgEjAS6AGMF5EetXa7GSgwxnQCngIe93WhNfL37WNAyQK2pF8EkbEN9TRKKXVqWvaDm7+CsS/bGRxfGQ6f3GaHPDYQb1roA4CNxpjNxphyYCpQu2k8Bnjduf4hMEykYVaZ2PTN+zSRUmIztbtFKeXnRKDPVXY2x3Pvg9Uf2W6YVR82yNN5E+itgB0et3Oc++rcxxhTCRQCKbUfSEQmisgSEVmSl5d3UgVHxiexLO5s2vcfdlI/r5RSjS46AYY/ArcvhPYXQEqnBnmaiAZ51GMwxrwEvASQmZl5UoM1+w0fD8N13nOlVABK6WiHOTYQb1roO4E2HrdbO/fVuY+IRABJQL4vClRKKeUdbwJ9MdBZRNqLSBQwDphWa59pwPXO9SuA2cYE0KzwSikVBOrtcjHGVIrIncAMIBx41RizRkQeBZYYY6YBrwBvishGYB829JVSSjUir/rQjTHTgem17nvY43opoOfgK6WUiwLyTFGllFJH00BXSqkgoYGulFJBQgNdKaWChLg1ulBE8oBtJ/njqcBeH5YTCPSYQ4Mec2g4lWNuZ4ypc95w1wL9VIjIEmNMptt1NCY95tCgxxwaGuqYtctFKaWChAa6UkoFiUAN9JfcLsAFesyhQY85NDTIMQdkH7pSSqmjBWoLXSmlVC0a6EopFST8OtD9bXHqxuDFMd8vImtFZKWIzBKRdm7U6Uv1HbPHfpeLiBGRgB/i5s0xi8hVzu96jYg03KoIjcSL93ZbEZkjIsuc9/coN+r0FRF5VURyRWT1MbaLiDzjvB4rRaT/KT+pMcYvL9ipejcBHYAoYAXQo9Y+twMvONfHAe+5XXcjHPMQIM65/ptQOGZnvwRgPrAQyHS77kb4PXcGlgHNnNvpbtfdCMf8EvAb53oPYKvbdZ/iMZ8P9AdWH2P7KOBLQIBBwKJTfU5/bqH71eLUjaTeYzbGzDHGFDs3F2JXkApk3vyeAf4CPA6UNmZxDcSbY74VmGSMKQAwxuQ2co2+5s0xGyDRuZ4E/NSI9fmcMWY+dn2IYxkDvGGshUBTETntVJ7TnwPdZ4tTBxBvjtnTzdi/8IGs3mN2Poq2McZ80ZiFNSBvfs9dgC4i8q2ILBSREY1WXcPw5pgfASaISA52/YW7Gqc015zo//d6Neoi0cp3RGQCkAlc4HYtDUlEwoAngRtcLqWxRWC7XQZjP4XNF5Hexpj9rlbVsMYDU4wxT4jIWdhV0HoZY6rdLixQ+HMLPRQXp/bmmBGR4cBDwGhjTFkj1dZQ6jvmBKAXMFdEtmL7GqcF+Bej3vyec4BpxpgKY8wWYD024AOVN8d8M/A+gDHmeyAGO4lVsPLq//uJ8OdAD8XFqes9ZhHpB7yIDfNA71eFeo7ZGFNojEk1xmQYYzKw3xuMNsYscadcn/Dmvf0ptnWOiKRiu2A2N2aRPubNMW8HhgGISHdsoOc1apWNaxrwK2e0yyCg0Biz65Qe0e1vguv5lngUtmWyCXjIue9R7H9osL/wD4CNwA9AB7drboRj/hrYAyx3LtPcrrmhj7nWvnMJ8FEuXv6eBdvVtBZYBYxzu+ZGOOYewLfYETDLgYvcrvkUj/ddYBdQgf3EdTNwG3Cbx+94kvN6rPLF+1pP/VdKqSDhz10uSimlToAGulJKBQkNdKWUChIa6EopFSQ00JVSKkhooCulVJDQQFdKqSDx/wG9v0SqTrAfVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Normalized spline & difference curve\");\n",
    "plt.plot(x_normalized, y_normalized);\n",
    "plt.plot(x_difference, y_difference);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: Identify local maxima  \n",
    "of the difference curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import argrelextrema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# local maxima for knees\n",
    "maxima_indices = argrelextrema(y_difference, np.greater)[0]\n",
    "x_difference_maxima = x_difference[maxima_indices]\n",
    "y_difference_maxima = y_difference[maxima_indices]\n",
    "\n",
    "# local minima\n",
    "minima_indices = argrelextrema(y_difference, np.less)[0]\n",
    "x_difference_minima = x_difference[minima_indices]\n",
    "y_difference_minima = y_difference[minima_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8dcn+04gC0sgBAKI4MISERBXUBEXrCtYd1zq2n7tZq21rbX9dXXHKm5YbUWr1VJBtG6VBBEQEGQzYQ1rMiGBhOzJ+f1xbmAICRlgkjsz+Twfj3kwc+/N3M9lknduzj3nXDHGoJRSKviFuV2AUkop/9BAV0qpEKGBrpRSIUIDXSmlQoQGulJKhQgNdKWUChEa6AFGRDaJyIR2fP+zRGRre73/EdTxrIj8ogP2kykiFSIS7qf3MyIywHl+0DGIyB0issvZX4qInCYi+c7rS/2xf6UOJ8LtAlTnZIz5XgftZwuQ0E7vvf8YRCQSeBQYbYz52ln2MPC0MeaJ9ti/Us3pGbpS/tEdiAFWeS3r2+y1z0QkpE62Qu14ApUGegATkWgReVxEtjuPx0Uk2mv9ZBFZLiJ7RWS9iEx0lt8kImtEpFxENojI7UewTyMidzpNBeUi8hsRyRaRBc5+3hSRKGfbriLynogUi0ip87y3s66biGwVkYud1wkiUiAi1zuvZ4rII87zs5xtfyIiRSKyQ0QuFZFJIvKtiOwWkQe8ahwlIl+ISJmz7dNNNbVwPFnOMUU4rz9zjinPOb4PRST1MP8fP3b2sV1Ebm62bqaIPCIig4B1zuIyEflERNYD/YH/OE0u0SLSRURedN5vm/O14c573ejU9JiIlAC/cpbf7HyWpSLygYj0bfZZfc/5rMpEZLqIiNf6W72+D1aLyAhneS8Redv53DaKyL2HOf5YEfmLiGwWkT0ikussO6TpTryaC0XkVyLyloi8JiJ7gQdEpEpEunltP1xEPGL/ujnssSofGWP0EUAPYBMwwXn+MLAQSAfSgAXAb5x1o4A9wLnYX8wZwGBn3YVANiDAmUAlMMJZdxaw9TD7N8C/gSRgKFADfIwNpy7AauAGZ9sU4HIgDkgE/gm86/Ve5wE7nfqfB97yWjcTeMSrpnrgISASuBUoBv7hvO9QoAro52w/EhiNbTLMAtYAP2jleLKcY4pwXn8GrAcGAbHO69+38rUTgV3ACUC8U48BBrRwDAftp/ln6bx+B3jOea90YBFwu7PuRuf/4B7nuGKByUABcLyz7EFgQbPP6j0gGch0/s8mOuuuBLYBpzjfBwOwfzGEAV85/9dRzue6ATi/lf+D6c7/UQYQDowFomnh+4iDv3d/BdQBlzr7jAU+AW712v5PwLPO88Meqz58zA+3C9BHsw/k4B+K9cAkr3XnA5uc588Bj/n4nu8C33eeH/KD2GxbA5zm9for4Kder/8CPN7K1w4DSpstewpY6YRLitdy7zA8CxvY4c7rRKeOU5vVcWkr+/0B8E4r67I4NNAf9Fp/JzCvla99Ca+wx/4SOKpAxzbJ1ACxXuunAp86z28EtjTb//vANK/XYdhfzn29PqtxXuvfBO53nn/Q9Jk3e89TW9jPz4CXW9g2zPlcTm5h3SHfRxwa6J83W38L8InzXIBC4AxfjlUfvj20ySWw9QI2e73e7CwD6IMN/EOIyAUistBpqigDJgGtNiu0YJfX86oWXic4+4kTkeecP8f3Ap8DyXJwj5IZ2DPcmcaYksPss8QY0+C1j5bqaNrvIKd5Z6ez398d4fHt9HpeSesXTXthQ6fJ5la280Vf7F8fO5zmkTLsL+V0r20KW/iaJ7y2340NwgyvbVo7lta+P/oCvZre03nfB7C/cJpLxV4XaPH7zAfNj+dtYIyI9ATOABqB+V51tXWsqg0a6IFtO/YbvUmmswzsD0t28y8Q28b+NvBnoLsxJhmYi/3h8LcfAsdhz6STsD+kNO3LCfYZwN+AO8Xp7ucHfwXWAgOd/T5A+xzfDmwwNsk8hvcqxJ6hpxpjkp1HkjFmqNc2zac+LcQ2ySR7PWKNMQt83N8h3x/O8o3N3jPRGDOphW09QHUr77MP29QG7P+s05ptc9DxGGNKgQ+Bq4FrgFnGOR3n2I5VOTTQA9vrwIMikuZcuHsIeM1Z9yJwk4iMF5EwEckQkcHYdtFobHtqvYhcgG3Lbg+J2DPnMudi1y+brX8A+0N9M7a99G/in/7gicBeoMI55jv88J4teRO4UUSGiEgchx6fz4wxO7Bh9hcRSXI+s2wROfMwX/Ys8DMRGQrgXFS90sddvgD8SERGijXAuci4CCgXkZ86FzfDReQEETmlhZobsc1OjzoXUsNFZIxz0vAtECMiFzoXNR/Eft+15R/A9cAVznN/HKtyaKAHtkeAJcAKbDv0UmcZxphFwE3AY9iLo//DtjeWA/diw6gUeyY0u53qexx7scuDvXg7r2mFiIwE7gOud5pS/oAN9/v9sN8fYY+rHHux9Q0/vOchjDHvY4/xE+wFu0+O8S2vx/7CXY39bN4Ceh5m/+9g/99mOU1L3wAX+LIjY8w/gd9iQ7Mcex2lm/NZXIS93rER+9m9gL3g3ZIfYb/3FmObQf4AhBlj9mCvP7yAvT6yD/BlwNpsYCCw0zj99Y/1WNUBcuAvHqWUUsFMz9CVUipEaKArpVSI0EBXSqkQoYGulFIhwrUJc1JTU01WVpZbu1dKqaD01VdfeYwxzfv8Ay4GelZWFkuWLHFr90opFZREpNURy9rkopRSIUIDXSmlQoQGulJKhQgNdKWUChEa6EopFSLaDHQReUnsbcG+aWW9iMiTYm8vtqLpNldKKaU6li9n6DOxt+JqzQXY2dMGArdh56pWSinVwdrsh26M+VxEsg6zyWTgb85E9QtFJFlEejrzPyulVNAyxlBRU09ZZZ19VNVSWlnHnkr77wkZSZwzuKWbPbnDHwOLMjj4VlNbnWWHBLqI3IY9iycz81hu/qKUUr4zxlBV10BZZR2llbXsqayj1AloG9Y2oJuel1U5/1bWUd/Y+hTjN47NCrlA95kxZgb2lmTk5OToROxKqSNmjKG8ph5PeQ2eilp276txgtoJ6H3eZ9I2wMuq6qitb2z1PWMjw+kaF0mXuCi6xkUyqHsCyXFRJMdG0jUuii5x9t/kuEi7XWwUXWIjiYoIrH4l/gj0bRx838XezjKllPKJMYa91fV4KmrwlNdQ7Pzrqai1yypqKK6odZbVUNNKOEdFhNE1LpLkWBu+WalxDI9LPhDIsZE2qL0CuktsJDGR/rgzovv8EeizgbtFZBZwKrBH28+VUsYY9lbVU1xRTXH5gWC2oV1rQ7spuPfVtngGHR4mdIuPIjUhmtSEKLLT4klLiLavE+3ybvFRdI2zj5jIMETa437hwaHNQBeR14GzgFQR2Yq9UW4kgDHmWewd5Sdh77lYib3PpVIqRDWdTW8vq9r/2LGnmuLypsC24V1SUUttQ8shneKEdFpiNAPTE0lNjDoQ1M7y1AQb0mFhnTegj5QvvVymtrHeAHf5rSKllKtq6xvZtbeabV6Bva2s2gnuKraXVVNRU3/Q10SECSkJUU4QR3Ncj8T9Z9VpidE2rJ11ybGRGtLtxLXpc5VSHc8YQ2llnRPSBwJ7+57q/c+Lymtofu/4lPgoeiXH0i81ntMGpJKRHEvPLrH0So4hIzmW1IRoDekAoIGuVAiprmtghxPO28qq2OGcWW/fcyDAq+sObgaJjgizAZ0cwxkD0+iVHEtGciy9km1g9+wSS2xUaFw0DHUa6EoFmYZGw/ayKgqKK1hfVMH64n2sL65gQ/E+PBU1h2yflhhNr+RYBvdI5Jzj0p2gbgrtGLrFR3XqC4mhRANdqQBVWVvPBiesm0J7fVEFGz37Duq2lxwXyYC0BM4ZnEafrnH0cs62M5Jj6dElhugIPbvuLDTQlXKRMYbi8hp7tl28zznjtmfb28qq9m8XJtCnWxzZaQmcPjCV7LQEstMTyE5LoFt8lItHoAKJBrpSHaC2vpEtu/dRUNR0xm0DfENRBeVePUbiosLJTkvglKyuTEnrsz+0+6bEhczgF9V+NNCV8qPy6jryiw5u215fXMGWksqD5gTpnhRNdloClw7PIDstngHpiWSnx9MjKUbbs9VR00BX6ig1NhryiypYtqWUZVvKWLqllILiiv1d/iLDhayUeAalJ3LBCT1sM0laAv3T4kmMiXS3eBWSNNCV8lHpvlqWF5axbEspS7eU8XVh2f7mkuS4SIb3Seaik3oxtFcS2ekJ9OkaS0R4YE3epEKbBrpSLahvaGTdrvL9Z97Lt5SxwbMPsBcoB/dI4pJhvRiR2ZXhmcn0S43XphLlOg10pQBPRc3+8F62pZQVW/dQWdsAQGpCFMMzu3JFTm+G9+nKSb27EB+tPzoq8Oh3pep0ausbWbtzL0s3l7Ks0IZ44W7bRTAiTBjaK4mrcvowPDOZEZld6d01Vs++VVDQQFchb9fe6gPhvbmUldv27B+Y0z0pmhGZXbl+dBbDM5M5IaOLdg9UQUsDXYUUYwzrdpWTm+9h2RZ7AXP7nmrA3vzgxIwuXDe6L8MzuzKibzI9u8S6XLFS/qOBroKeMYa1O8uZu3IHc1bs2H/xsnfXWEZmdePWzGSGZ3bl+J6JOgxehTQNdBWUWgrxMIHR/VOYdno/Jhzfne5JMW6XqVSH0kBXQaMpxOes2MHclYeG+PlDe5CaEO12mUq5RgNdBTRjDGt22DNx7xAfk60hrlRzGugq4HiH+JyVO9ioIa6UTzTQVUA4XIjfoiGulE800JVrmkJ8zsrtzF25U0NcqWOkga46lDGG1Tv2Om3iB4f4raf35/yh3UnREFfqqGigq3bnHeJzVuxgU0mlhrhS7UADXbWL1kJ8bHYqt52RrSGuVDvQQFd+t3bnXn47Zw3z8z0a4kp1IA105TdF5dU89t9veWNxIYkxkfx80vFcNiJDQ1ypDqKBro5ZdV0DL+Zu5JlPC6ipb+TGsf24d/wAkuP0bvRKdSQNdHXUjDHM/no7f5y3jm1lVZw3pDv3XzCY/mkJbpemVKekga6Oylebd/Ob99awvLCMob2S+POVJzMmO8XtspTq1DTQ1REp3F3J7+etZc6KHaQnRvOnK07i8hG9CQvTO/oo5TYNdOWTvdV1TP+0gJdzNxEWBt8fP5Dbz+xPXJR+CykVKHz6aRSRicATQDjwgjHm983WZwKvAMnONvcbY+b6uVblgvqGRl5fXMhj//2W3ftquXxEb358/nH06KJzjSsVaNoMdBEJB6YD5wJbgcUiMtsYs9prsweBN40xfxWRIcBcIKsd6lUd6NN1Rfxuzhryiyo4tV83HrxwCCf27uJ2WUqpVvhyhj4KKDDGbAAQkVnAZMA70A2Q5DzvAmz3Z5GqY63bWc5v567h82+LyUqJ47nrRnLekO6I3vleqYDmS6BnAIVer7cCpzbb5lfAhyJyDxAPTGjpjUTkNuA2gMzMzCOtVbWz4vIaHvvoW2Yt2kJCdAS/uGgI143uS1REmNulKaV84K8rWlOBmcaYv4jIGOBVETnBGNPovZExZgYwAyAnJ8f4ad/qGFXXNfBS3kae+XQ91XUN3DA2i3vPGUjXeB0YpFQw8SXQtwF9vF73dpZ5mwZMBDDGfCEiMUAqUOSPIlX7MMbwnxU7+MP7a9lWVsWE47vzwCQdGKRUsPIl0BcDA0WkHzbIpwDXNNtmCzAemCkixwMxQLE/C1X+9dXmUh6Zs5plW8oY0jOJP11xEmMHpLpdllLqGLQZ6MaYehG5G/gA2yXxJWPMKhF5GFhijJkN/BB4XkT+D3uB9EZjjDapBKDC3ZX8Yd5a3nMGBv3RGRgUrgODlAp6PrWhO33K5zZb9pDX89XAaf4tTflTeXUd0z9dz0t5GwkTuHf8QG4/oz/x0TowSKlQoT/NIa6+oZE3lhTy6IffUrKvlstGZPDj84+jZ5dYt0tTSvmZBnoI27W3mhteWsTaneWM6teNly88npN6J7tdllKqnWigh6h9NfXcPHMxhbsr+et3RzDxhB46MEipEKeBHoLqGxq55/VlrN1Zzgs35HD2celul6SU6gA6BDDEGGP49X9W88naIn59yVANc6U6EQ30EPNi7kZeXbiZ28/oz7Wj+7pdjlKqA2mgh5D3V+7gt3PXMOnEHvx04mC3y1FKdTAN9BCxbEspP3hjOcP6JPPoVcP0DkJKdUIa6CFgS0klt7yyhO5JMbxwfQ4xkeFul6SUcoEGepArq6zlxpmLqG80vHzTKaQkRLtdklLKJRroQaymvoHbX/2KrburmHHdSLJ1lkSlOjXthx6kjDHc//ZKvty4myemDOPU/ilul6SUcpmeoQepxz7K551l2/jReYOYPCzD7XKUUgFAAz0I/XNJIU9+nM9VOb256+wBbpejlAoQGuhBJq/Aw8/+tZJxA1L57XdO1PlZlFL7aaAHkW93lfO9176if1o8z1w7gshw/fiUUgdoIgSJovJqbnp5MTGR4bx80yiSYiLdLkkpFWA00INAZW09t7yyhN37annphlPISNabUyilDqWBHuAaGg33vr6cb7bt4ampwzmxdxe3S1JKBSgN9AD3yJzVfLRmF7+8eCgThnR3uxylVADTQA9gL+dt5OW8TUwb148bxma5XY5SKsBpoAeoD1ft5OH3VnP+0O48MOl4t8tRSgUBDfQA9HVhGffOWsZJGV14/OrhhOtUuEopH2igB5jC3ZVMe2UJqQnRvHDDKcRG6VS4Sinf6ORcAWRPVR03z1xMbX0Ds247lbREnQpXKeU7DfQAUVvfyB2vfcWmkn28cvMoBqQnul2SUirIaKAHAGMMD7yzkgXrS/jLlSczNjvV7ZKUUkFI29ADwFOfFPDWV1v5wYSBXD6yt9vlKKWClAa6y95ZtpVH//stl43I4PvjB7pdjlIqiGmgu2jhhhJ+8tYKxvRP4feXnaRT4SqljolPgS4iE0VknYgUiMj9rWxzlYisFpFVIvIP/5YZegqKKrjtb0vomxLPs9eOJCpCf7cqpY5NmxdFRSQcmA6cC2wFFovIbGPMaq9tBgI/A04zxpSKSHp7FRwKPBU13DRzEVERYbx84yl0idOpcJVSx86X08JRQIExZoMxphaYBUxuts2twHRjTCmAMabIv2WGjqraBm55ZQnF5TW8cMMp9OkW53ZJSqkQ4UugZwCFXq+3Osu8DQIGiUieiCwUkYktvZGI3CYiS0RkSXFx8dFVHMQaGw3/98Zyvt5axhNThjOsT7LbJSmlQoi/Gm4jgIHAWcBU4HkROSStjDEzjDE5xpictLQ0P+06ePy/99cwb9VOHrxwCOcP7eF2OUqpEONLoG8D+ni97u0s87YVmG2MqTPGbAS+xQa8crz6xSaen7+RG8dmcfNpWW6Xo5QKQb4E+mJgoIj0E5EoYAowu9k272LPzhGRVGwTzAY/1hnUPlm7i1/OXsWE49P5xUVDtHuiUqpdtBnoxph64G7gA2AN8KYxZpWIPCwilzibfQCUiMhq4FPgx8aYkvYqOpg0Nhruf3slg3sk8eRUnQpXKdV+fJrLxRgzF5jbbNlDXs8NcJ/zUF7W7iynqLyGn0wcTFyUTp2jlGo/OpqlneUVeAAYN0An3FJKtS8N9HY2v8DDgPQEenSJcbsUpVSI00BvR9V1DSzaWKJn50qpDqGB3o6Wbi6luq6R0wdqoCul2p8GejuaX+AhIkw4tX+K26UopToBDfR2lFfgYXhmMgnR2rtFKdX+gjJpzjrrLLdLaFNDRAyFI+8meWseZ816oMP3Hx/ewKSeJUzovpt99eEUVsZQWBXNlspotlbGsLM6ika0T7xSbvjss8/a5X2DMtCDQXVSJogQs2dzh+63V0wNl/UuZlKPEuIiGlm9N47osEbOTi8lKbJh/3a1jcK2qmgKK6P3h33T8731+m2hVDAKyp/c9vrt5k8/+9cK3luxg9x/v0ZEeDu3bBkDmxfAwmdg7RwIC4cTroTRdzCk1/AD21SWgCcfSgqIKsmnn6eAfiX5sHsjNNYdeL/YbpA6EFIG2EfqQEgZCN36QUR0+x6LUuqoBWWgB4PcAg9j+qe0b5jX18Kqf9kg3/E1xHaF0++DU26FpJ4HbysC8an20XfMwesa6qFsM5QUOIGfD54CKPgYlv/d6z3CIDnThntT4Df9m9jT7kMp5RoN9HawuWQfhburuPX0/u2zg30l8NVLsOgFqNgJqcfBRY/DSVdD1FHcMCM8AlKy7WPQ+Qevq95rg77p0RT4m/OgrvLAdlEJznt4hX36EEg/XoNeqQ6igd4O5ue303D/4nX2bPzrWVBfDdnnwOTp9t+wdvpLICYJMkbYh7fGRijf7hXyzr9bF8E3bwPGbpc+BIZfa3/ZxGt/fKXakwZ6O8jN95CRHEu/1PhjfzNjYP3HsPCvUPARhEfDyVfD6Dvt2a9bwsKgS2/76H/WwevqqmH3BtjyBSz/B3zwAPz3IRg00Yb7gHPtXwVKKb/Snyo/a2g0LFjv4YITeh7bvOd1VbDiDRvkxWshoTuc/SDk3BT4Z7qRMdB9iH2cMg2K1sLy1+xfFmvfs8dy8hQYdi2kDXK7WqVChga6n63ctoe91fWcdrTD/ct3wuIXYMlLtldKjxPh0mfhhMuCt4dJ+mA47xEY/0vI/xCW/R0WPA15T0DvUTD8uzD0Mtu8o5Q6ahrofpabb29+fVr2EQ733/E1fPGMbX9urIfjJsGYO6HvaaFzUTE8EgZfaB8VRfYvkGWvwX++D+/fD0MvhWHftcfcXtcElAphGuh+Nj/fw9BeSaQk+HA23dgA386zQb45FyLjIedmOPV222MklCWkw9h7YMzdsO0rG+zfvA1fvw5ds2ywnzwVkvu0+VZKKUsD3Y/21dSzdEspN4/rd/gNayps/+6Ff4XSjdClj22SGH4dxCZ3TLGBQgR659jH+b+DNf+x7e2f/hY+/R1kn23DffBFtm1eKdUqDXQ/WrRpN3UNpvXuimVb4MvnYOmrULPHth9P+CUMvlh7fYDtQ3/y1fZRusn2kFn+D3h7GsR0gROvtL1keg4LnWYopfxIU8SPcvM9REWEcUpWt4NX7PwGPv+TPfsEGDIZxtxlz0pVy7pmwdkPwJn3w8b/2b9olr5qLxh3P8GetZ90VeD3+FGqA2mg+1FuvodRWd2IiQw/sLChDv422c6VMvZuOyxf24V9FxZmm12yz4ZJpbadfdnf4YOf2b7tx020TVXZ4/WvHNXp6U+AnxTtrWbdrnK+MyLj4BUbP4dKD0x5HQZPcqe4UBHbFU65xT52rbLBvuIN+5dPQg/bt334tXbqAaU6Ie0b5id561sZ7r/6XYhKtMPzlf90HwoTfwf3rYGrX4New2DBU/B0Drx4nh3E1FDvdpVKdSgNdD+Zn++hW3wUQ3p6DY5pqLNnj8ddoD002ktEFBx/MVzzBty3Gib8Gip3wzu3w9Mjbbt7Q13b76NUCNBA9wNjDLn5HsZmpxAW5tX7YuPnUFVqB8yo9pfYA8b9AO5eDFNnQUwyzL4bnhoBX8200w0rFcI00P0gv6iCovIaTm8+3H9/c8t4dwrrrETsX0W3fQbX/BPi0+xo1KdGwOIXob7G7QqVahca6H6wf7rcgWkHFjbUwZr3bC8MbW5xhwgMOg9u+RiufdvehGPOffDkcFj0vJ0VUqkQooHuB3kFHvqlxpORHHtg4ab5ULUbhmhzi+tEYMAEmPYhXPeOHZk790fw5DBY+Kyd2VKpEKCBfoxq6xtZuKHk0N4tq961d/EZoM0tAUPE9ja6eR5cPxu69Yd5P4UnToYvpkNtZdvvoVQA00A/Rsu2lFJZ28A47/bzhno77/egiRAZ2/oXK3eIQP8z4aa5cOMcSB1kb8LxxEmQ9yTU7nO7QqWOigb6Mcot8BAeJozxni5303w7l/nQ77hXmPJN1ji48T246X3bt/2/v4DHT4Tcx+wkakoFEZ8CXUQmisg6ESkQkfsPs93lImJEpNNMUpJb4OHk3l1Iiok8sHC1NrcEnb5j4fp/w80f2sm/PvqVDfbP/2xvlK1UEGgz0EUkHJgOXAAMAaaKyJAWtksEvg986e8iA9Weqjq+Liw7uP28od4OJtLmluCUeSpc9y/bM6Z3DnzyGxvs//sjVO9xuzqlDsuXM/RRQIExZoMxphaYBUxuYbvfAH8AOk1fsC/Wl9BomnVX3N/cor1bglrvHPjuP+HWTyFzjJ2f/bET4dP/ZweLKRWAfAn0DKDQ6/VWZ9l+IjIC6GOMmXO4NxKR20RkiYgsKS4uPuJiA01uQTHxUeEMz/S6KcX+5pYJ7hWm/CdjBFwzC27/HPqdDv/7PTx+EnzyiJ1iQKkAcswXRUUkDHgU+GFb2xpjZhhjcowxOWlpaW1tHvDyCkoY3T+FyHDnv3F/c8v52twSanqeDFP+Dt/Lhf5n2fntHz8RPvo17CtxuzqlAN8CfRvgPYF3b2dZk0TgBOAzEdkEjAZmh/qF0a2llWz07OM07/bzzbm2uUUHE4WuHifC1a/CHV/AwHNtb5jHT7Rzs+/zuF2d6uR8CfTFwEAR6SciUcAUYHbTSmPMHmNMqjEmyxiTBSwELjHGLGmXigNErjPc/6D5W1a9a2/0PPBcl6pSHab7ELhyJty50M4bk/ekDfYPfg4Vwd+cqIJTm4FujKkH7gY+ANYAbxpjVonIwyJySXsXGKjmF3jonhTNgPQEu0CbWzqn9MFwxYtw1yI7je/CZ+xcMfP/olMKqA7nUxu6MWauMWaQMSbbGPNbZ9lDxpjZLWx7VqifnTc2GhYUeBg3IA1pulnx5jx7ZyLt3dI5pQ2Cy2bAnV/ai6cfPwxP5cDXb0Bjo9vVqU5CR4oehdU79lJaWXdwc8tqp7llgDa3dGppg2Dq63DDexCfAu/cBi+cA5vy3K5MdQIa6EehabrcsQOc4f7ezS1RcS5WpgJGv9Ph1s/gO89BRRHMnASzvgueArcrUyFMA/0o5BYUM7hHIumJzjznm/NgX7E2t6iDhYXZG1ffvQTOeRA2fAbPnArv/1T7sKt2oYF+hKrrGli8qfTg4f6r34XIOG1uUS2LivqttD8AABF0SURBVIMzfgz3LoPh18KiGfDEMNszRu+epPxIA/0ILd60m9r6xgPT5TY2aHOL8k1COlz8BNyxAPqMsjM7Pn0KfPMvMMbt6lQI0EA/Qrn5HqLCwxjVr5td0NTcooOJlK/Sj4dr37J3T4pKgLdughfPg8JFblemgpwG+hGan+9hRN9k4qIi7IJVTnPLwPPcLUwFn+xz4Hvz4ZKnoGwzvHgu/PNG2L3R7cpUkNJAPwKeihpW79jL6U2zKzY1tww8T5tb1NEJC4cR18M9S+HMn8K6eTB9lB1xqrM6qiOkgX4E8gpsd8X9F0Q3L4B9Rdq7RR276AQ4+wG4dymceJW9x+mTw+1NrOtr3a5OBQkN9COQV+ChS2wkJ2R0sQtWa3OL8rOkXnDpdDtdb4+T7E2snxlt/xLUC6eqDRroPjLGkJvvYWx2CuFhYptbVs92mlvi3S5PhZqeJ9lb4l3zJoRFwBvXwswLYdtStytTAUwD3UcbPPvYvqf6QHdFbW5R7U3Edoe9YwFc+CgUr4Pnz4a3b4Wywra/XnU6Gug+2j9d7gDngujqdyEiVptbVPsLj4BTptmBSePug9X/hqdz7M019AbWyosGuo9yCzxkdosjMyXuQHPLIG1uUR0oJgkm/BLu+QqOvwRyH7UXThe/YOcTUp2eBroP6hsaWbi+5MDdibZ8YZtbdDCRckNyH7j8ebj1E0gdBHN+CH8dC99+oBdOOzkNdB98vbWM8pr6A9PlrnKaWwad725hqnPLGAk3zYWr/w6N9fCPq+DvV0DJercrUy7RQPfB/HwPIjA2O8UZTDTb3mZOm1uU20Tg+IvsrfDO/x1s+dJ2c/z4Yajd53Z1qoNpoPsgr8DDSRldSI6Lgi0LoWKX9m5RgSUiCsbcBfcsgaHfsbfAe3qUvYCqzTCdhgZ6Gypq6lm2pexA+/nqdyEiBgZqc4sKQIk97K3wbnofYrrAm9fDq98BT77blakOoIHehoXrS6hvNLb/ufdgougEt0tTqnV9x9rRphP/ANu+gmfGwH9/CTUVblem2pEGehtyCzzERoYzsm9Xp7llpza3qOAQHgGjv2e7OZ54JeQ9bif+0vnXQ5YGehvm5xczql83oiPCtblFBaeEdPjOX+HmDyGum51//W+T7chTFVI00A9jx54q1hfvs90VGxud5pZztblFBafMU+G2/8GkP8OO5bbv+ocPQk2525UpP9FAP4ym4f6nDUiFQqe5RQcTqWAWFg6jbrXzr588BRY8ZW+Dt/ItbYYJARroh5Fb4CE1IZrBPRKdwUQxMGii22UpdeziU2HydJj2kW2SeXsavHIx7FrtdmXqGGigt6Kx0ZBX4GHcgBTEGDuYaMAEbW5RoaXPKXDrp3DRY7BzJTw7DuY9oJN+BSkN9Fas21WOp6KWcQPToPBLKN9hB2woFWrCwiHnZtsMM+I6WPiMnc3x6ze0GSbIaKC3oqn9fNyAVNu7JTxa525RoS0+BS5+Am79GJIy4J3b4OVJsPMbtytTPtJAb8X8Ag8D0hPokRhlh08PPBeiE90uS6n2lzESbvkYLn4SitfCc2fA+z+FqjK3K1Nt0EBvQXVdA4s2ltizc21uUZ1RWBiMvMEOShp5I3z5nG2GWf4P24VXBSSfAl1EJorIOhEpEJH7W1h/n4isFpEVIvKxiPT1f6kdZ+mWUqrrGm3/c21uUZ1ZXDe46FG47TPomgXv3gEvT4QdK1wuTLWkzUAXkXBgOnABMASYKiJDmm22DMgxxpwEvAX80d+FdqTcfA8RYcKp/bpqc4tSAL2G2ZGmk6fb+dZnnAlzfgRVpW5Xprz4coY+CigwxmwwxtQCs4DJ3hsYYz41xlQ6LxcCvf1bZsfKLfAwPDOZhKKltrlFBxMpZZthhl9rp+g95RZY8iI8NRKWvqrNMAHCl0DPALxvMb7VWdaaacD7La0QkdtEZImILCkuLva9yg5Uuq+Wldv2MG5Amh1MFB4Nx+lgIqX2i+0Kk/5kpxFIGQiz74aXztNmmADg14uiInItkAP8qaX1xpgZxpgcY0xOWlqaP3ftNwvWl2AMjBvQzTa3DJigzS1KtaTnSXDzPLj0Wdi9EWacZeeG0TslucaXQN8G9PF63dtZdhARmQD8HLjEGFPjn/I6Xm6Bh8ToCE7mWyjfrlPlKnU4IjBsKty92DbHLHgKpp8K6+a5XVmn5EugLwYGikg/EYkCpgCzvTcQkeHAc9gwL/J/mR0nt6CY0dkpRKyd7fRu0eYWpdoU1w0ueRJummfvtfv61fDGdbB3u9uVdSptBroxph64G/gAWAO8aYxZJSIPi8glzmZ/AhKAf4rIchGZ3crbBbTNJfso3F3F6fubW8ZDTJLbZSkVPPqOgdvnw/iHIP9De1/TL5+zd/tS7S7Cl42MMXOBuc2WPeT1fIKf63LFfGe4//jELbB3G0z4lav1KBWUIqLg9B/awXhzfgjv/wS+ft1OK9DzZLerC2k6UtRLXoGHjORYem2dp80tSh2rbv3h2n/B5S/Cnq32oukHP9f7mrYjDXRHQ6NhwfoSxmV3RdbM1uYWpfxBBE68wl40HXEDfPG0vWi6dm7bX6uOmAa6Y+W2PeypquOilO22uUUHEynlP7Fd4eLH7WjTmCSYNRVmfRf2HNJhTh0DDXRHbr4d6JSz738QHqWDiZRqD5mnwu2f2+tTBR/D9FGw8Fm9aOonGuiO+fkeTuiZQGz+HMgeDzFd3C5JqdAUHgnj/g/u/AL6nArzfgrPnwPbl7ldWdDTQAcqa+tZuqWUK3sWwd6tOphIqY7QrR9c+zZc8ZKdM+n5c+D9+6Gm3O3KgpYGOvDlxt3UNRjGNy5wmlsucLskpToHETjhcrhrEYy8Cb581l40XfOe25UFJQ107HS5URFCxvYPIfscbW5RqqPFJtt516d9CDHJ8MZ34fVrbHdH5TMNdGygT+1ZhOzdqr1blHJTn1Fw+/9gwq9h/Sd2pOkX06Gh3u3KgkKnD/Si8mrW7SrnspjF2tyiVCAIj4RxP4C7FkLWafDBA/D82bBtqduVBbxOH+h5BR7AcHzpp7a5JTbZ7ZKUUmBveXfNm3DlK1BRBC+Mh7k/geq9blcWsDp9oM/P93BG3GaiKnQwkVIBR8T2Ort7EeRMg0UzbN/11bPBGLerCzidOtCNMeTme7ihy3IIi9TmFqUCVUwXuPDPcMtHEJcKb14Hr0+Bsi1uVxZQOnWgFxRVUFRezejqXG1uUSoY9M6B2z6D8x6BjZ/bLo4LntKLpo5OHejz8z2cLOuJr9I7EykVNMIjYOw9cNeXkHW6ve3dixNg50q3K3Ndpw703AIPUxOWOs0tk9wuRyl1JJIz4Zo37EjTskJ47kz46NdQV+12Za7ptIFeW9/Iwg0ezmchZJ+tzS1KBaOmkaZ3L4aTp0Duo/DsabAp1+3KXNFpA33ZllIG1OXTtW6n9m5RKtjFdYNLn4Hr3oGGOph5Ifzn+1C9x+3KOlSnDfS8Ag8XhS/EhEXCYG1uUSokZJ9jZ3Ecczcs/ZsdadqJ5oXptIE+P7+YS6MWI9ln28n3lVKhISoezv+t7eIYn2rnhXnzeijf5XZl7a5TBvqeqjrMtqWkNxZpc4tSoSpjpO3ieM4vYN08mH6KPWsP4QFJnTLQv1hfwsSwL2nU5halQlt4JJzxI7gjD9KHwux74JWLoWS925W1i04Z6Hn5xVwU/iX0P1ObW5TqDFIHwo1z4KLHYMfX8NexkPdEyA1I6pSBXvTtQnpLMWFDv+N2KUqpjhIWBjk32wFJ2ePhvw/ZWRx3fO12ZX7T6QJ9a2klw8v/R6NE6GAipTqjpF4w5e92FsfynTDjbPjvL6Guyu3KjlmnC/Tcb4uZFLaQqt7jbN9VpVTn0zSL411fwrCpkPe4bYbZ+LnblR2TThfom1ctIDOsmLjhl7tdilLKbXHdYPJ0uP7fYBrtBdPZ90BVmduVHZVOFeiNjYb0wnk0EI4MvsjtcpRSgaL/WXDHFzD2Xlj22oE514NMpwr01dv3ML4hj+L0MdrcopQ6WFQcnPcbuPVTSEi3c67P+i7s3eF2ZT7rVIG+Ztl829wyTJtblFKt6DXMhvqEX0HBR3bO9a9mBsWApE4V6NHfzqaecJKG6ehQpdRhhEfCuP+DOxZAz5PsRF9BMCDJp0AXkYkisk5ECkTk/hbWR4vIG876L0Uky9+FHqvq2nqG7f2MzUk52tyilPJNSjbc8B+4+EnYsQKeGQPzH7UzOgagNgNdRMKB6cAFwBBgqogMabbZNKDUGDMAeAz4g78LPVarl80nU4qoGzzZ7VKUUsFEBEbeYLs4DjoPPv61HZC0fZnblR3ClzP0UUCBMWaDMaYWmAU0T8XJwCvO87eA8SIi/ivz2FUv/xf1JozMsVe6XYpSKhgl9YSrX4OrXoWKYnj+HPj8z25XdRBfAj0DKPR6vdVZ1uI2xph6YA+Q0vyNROQ2EVkiIkuKi4uPruKjFNY1kyVpk4lLTu/Q/SqlQsyQS+zZ+vDroFs/t6s5SERH7swYMwOYAZCTk9Ohl4xHX/XjjtydUiqUxSbDJU+6XcUhfDlD3wb08Xrd21nW4jYiEgF0AUr8UaBSSinf+BLoi4GBItJPRKKAKUDzIVSzgRuc51cAnxgTBJ02lVIqhLTZ5GKMqReRu4EPgHDgJWPMKhF5GFhijJkNvAi8KiIFwG5s6CullOpAPrWhG2PmAnObLXvI63k1oN1HlFLKRZ1qpKhSSoUyDXSllAoRGuhKKRUiNNCVUipEiFu9C0WkGNjcwbtNBTwdvM+OEsrHBqF9fHpswcuN4+trjElraYVrge4GEVlijMlxu472EMrHBqF9fHpswSvQjk+bXJRSKkRooCulVIjobIE+w+0C2lEoHxuE9vHpsQWvgDq+TtWGrpRSoayznaErpVTI0kBXSqkQEZKBHgo3tW6ND8d2n4isFpEVIvKxiPR1o86j1dbxeW13uYgYEQmYLmNt8eXYROQq5/NbJSL/6Ogaj5YP35eZIvKpiCxzvjcnuVHn0RCRl0SkSES+aWW9iMiTzrGvEJERHV3jfsaYkHpgp/hdD/QHooCvgSHNtrkTeNZ5PgV4w+26/XhsZwNxzvM7guXYfD0+Z7tE4HNgIZDjdt1+/OwGAsuArs7rdLfr9uOxzQDucJ4PATa5XfcRHN8ZwAjgm1bWTwLeBwQYDXzpVq2heIYeEje1bkWbx2aM+dQYU+m8XIi9w1Sw8OWzA/gN8AeguiOLO0a+HNutwHRjTCmAMaaog2s8Wr4cmwGSnOddgO0dWN8xMcZ8jr3PQ2smA38z1kIgWUR6dkx1BwvFQPfbTa0DkC/H5m0a9swhWLR5fM6fs32MMXM6sjA/8OWzGwQMEpE8EVkoIhM7rLpj48ux/Qq4VkS2Yu+tcE/HlNYhjvTnst106E2iVccRkWuBHOBMt2vxFxEJAx4FbnS5lPYSgW12OQv7l9XnInKiMabM1ar8Yyow0xjzFxEZg73D2QnGmEa3CwsloXiGHso3tfbl2BCRCcDPgUuMMTUdVJs/tHV8icAJwGcisgnbXjk7SC6M+vLZbQVmG2PqjDEbgW+xAR/ofDm2acCbAMaYL4AY7MRWocCnn8uOEIqBHso3tW7z2ERkOPAcNsyDpQ22yWGPzxizxxiTaozJMsZkYa8RXGKMWeJOuUfEl+/Ld7Fn54hIKrYJZkNHFnmUfDm2LcB4ABE5HhvoxR1aZfuZDVzv9HYZDewxxuxwpRK3ryC301XpSdizm/XAz51lD2N/+MF+M/0TKAAWAf3drtmPx/YRsAtY7jxmu12zP4+v2bafESS9XHz87ATbpLQaWAlMcbtmPx7bECAP2wNmOXCe2zUfwbG9DuwA6rB/RU0Dvgd8z+tzm+4c+0o3vyd16L9SSoWIUGxyUUqpTkkDXSmlQoQGulJKhQgNdKWUChEa6EopFSI00JVSKkRooCulVIj4/1m3+Y9nxzRjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"local maxima in difference curve\");\n",
    "plt.plot(x_normalized, y_normalized);\n",
    "plt.plot(x_difference, y_difference);\n",
    "plt.hlines(y_difference_maxima, plt.xlim()[0], plt.xlim()[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 5: Calculate thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sensitivity parameter S\n",
    "# smaller values detect knees quicker\n",
    "S = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tmx = y_difference_maxima - (S * np.abs(np.diff(x_normalized).mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 6: knee finding algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If any difference value (xdj, ydj), where j > i, drops below the threshold y = T|mxi\n",
    "for (x|mxi, y|mxi) before the\n",
    "next local maximum in the difference curve is reached,\n",
    "Kneedle declares a knee at the x-value of the corresponding\n",
    "local maximum x = x|xi.  \n",
    "**If the difference values reach\n",
    "a local minimum and starts to increase before y = T|mxi\n",
    "is reached, we reset the threshold value to 0 and wait for\n",
    "another local maximum to be reached.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# artificially place a local max at the last item in the x_difference array\n",
    "maxima_indices = np.append(maxima_indices, len(x_difference) - 1)\n",
    "minima_indices = np.append(minima_indices, len(x_difference) - 1)\n",
    "\n",
    "# placeholder for which threshold region i is located in.\n",
    "maxima_threshold_index = 0\n",
    "minima_threshold_index = 0\n",
    "\n",
    "curve = 'concave'\n",
    "direction = 'increasing'\n",
    "\n",
    "all_knees = set()\n",
    "all_norm_knees = set()\n",
    "\n",
    "# traverse the difference curve\n",
    "for idx, i in enumerate(x_difference):\n",
    "\n",
    "    # reached the end of the curve\n",
    "    if i == 1.0:\n",
    "        break\n",
    "\n",
    "    # values in difference curve are at or after a local maximum\n",
    "    if idx >= maxima_indices[maxima_threshold_index]:\n",
    "        threshold = Tmx[maxima_threshold_index]\n",
    "        threshold_index = idx\n",
    "        maxima_threshold_index += 1\n",
    "\n",
    "    # values in difference curve are at or after a local minimum\n",
    "    if idx >= minima_indices[minima_threshold_index]:\n",
    "        threshold = 0.0\n",
    "        minima_threshold_index += 1\n",
    "\n",
    "    # Do not evaluate values in the difference curve before the first local maximum.\n",
    "    if idx < maxima_indices[0]:\n",
    "        continue\n",
    "\n",
    "    # evaluate the threshold\n",
    "    if y_difference[idx] < threshold:\n",
    "        if curve == 'convex':\n",
    "            if direction == 'decreasing':\n",
    "                knee = x[threshold_index]\n",
    "                all_knees.add(knee)\n",
    "                norm_knee = x_normalized[threshold_index]\n",
    "                all_norm_knees.add(norm_knee)\n",
    "            else:\n",
    "                knee = x[-(threshold_index + 1)]\n",
    "                all_knees.add(knee)\n",
    "                norm_knee = x_normalized[-(threshold_index + 1)]\n",
    "                all_norm_knees.add(norm_knee)\n",
    "\n",
    "        elif curve == 'concave':\n",
    "            if direction == 'decreasing':\n",
    "                knee = x[-(threshold_index + 1)]\n",
    "                all_knees.add(knee)\n",
    "                norm_knee = x_normalized[-(threshold_index + 1)]\n",
    "                all_norm_knees.add(norm_knee)\n",
    "            else:\n",
    "                knee = x[threshold_index]\n",
    "                all_knees.add(knee)\n",
    "                norm_knee = x_normalized[threshold_index]\n",
    "                all_norm_knees.add(norm_knee)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhV1dn38e+dmYQkQMIcIMzzJAEEqyCgxQmqUMUZq1JRbH3U99Ha6tNq59ZZpOCE4gAItaJFLQSoA4OEeUgYZJCAZIQEEkKm9f6xDhAgkENyTnbOzv25rn2dk7N3su4Tkh8ra6+9thhjUEopFfiCnC5AKaWUb2igK6WUS2igK6WUS2igK6WUS2igK6WUS4Q41XB8fLxJTEx0qnmllApIa9asyTbGNK1sn2OBnpiYSEpKilPNK6VUQBKRvefap0MuSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhrobrd4sd2UUq7n2KX/qpb8/vf2cdQoZ+tQSvmd9tCVUsolNNCVUsolNNCVUsolqgx0EXlTRDJFZPM59ouIvCQiO0Vko4hc5PsylVJKVcWbk6IzgVeAd86x/yqgs2cbDEzzPKq6YPp0pytQStWSKgPdGPOliCSe55CxwDvGGAOsFJFGItLSGPODj2pUNdG1q9MVKBWQjDEUFpeRd6yEvGMl5HseK37cJ6ERo3o0d7rUk3wxbbE1sK/Cx+me184KdBGZBEwCaNu2rQ+aVlX65BP7eN11ztahlAPKyw1HjpeeDOOzQrnoxPPS04L6xHGl5ea8X/+uSxJdF+heM8bMAGYAJCUlnf87pXzj2Wftowa6CmDGGAqKy8g9WkxOwXFyC4rJKSgmt6CYQ4XFFQK79LTAPlJUwvkyOThIiIkIIbZBKLENQolpEEpC4wYnPz75ekToWa81jAghOEhq75vgBV8E+n6gTYWPEzyvKaVUpYwx5BeVkltQTG7BcXKOFp8W0qee2305BcUUl5ZX+rXCgoOIaRBKbAMbzHENw+jQNOqsQI6pGMiR9jEqLBiRuhXKNeGLQF8ATBGR2diToXk6fq5U/VJebsg7VlIhkI/b50dPhXROheA+VFhMSVnlXefIsGCaRIURFxVG04bhdG0eQ1zDsJOv2efhxEXZ1yJdFso1UWWgi8gHwHAgXkTSgf8DQgGMMf8AFgJXAzuBQuAufxWrlKp9hcWlHMwr4mB+ERn5RRzMO87BvGMczC/iYP5xMvKKyDp6nLJzjG1Eh4fQpKEN44TGkfRNaHTy4yZRJ4I6/ORrEaHBtfwO3cObWS43V7HfAA/4rCKlVK0oLzfkFhbbsD4tsE89/yGviCNFpWd9bnRECC1iImgRG0HnZvE0jwknLir8ZE/6REg3jgolPEQDurbo4lxuN2uW0xUoBxSVlJGZf9zTiy4iwxPSBys8Zh4pOmvYI0igaXQ4LWIb0D4+iiEd4mgeG2HD2xPgLWIjiAzT6KiL9F/F7dq0qfoYFXBKy8rZf/gYu7ML2JtTyJ6cAr7PKeRAnu1Z5xYUn/U5DUKDaRkbQfOYCAa1b0LzmAhaxNjwbuEJ7fiGYYQE64oggUoD3e3mzLGPN93kbB3qgpWUlZN+6Bh7cgrYm13AHk9w780pZF9u4WlzpCPDgmkXF0Wr2Aj6t21Ey5iIUz1rT4jHRIToyUOX00B3u2nT7KMGep1UUlbOvtxC9uYUenrbp4I7/dCx0040RoUFkxgfRY+WMVzduwXt4qJIjIsiMT6Spg3DNayVBrpS/lZcWs6+Q4XszSlgd/aJR9vT3n/49NBuGB5CYnwkvVrHcl2fVrSLi6R9fBTt4qKIbximoa3OSwNdKR/JOnKczQfy+C7z6MmhkT05Bew/dOy0qxWjw0NIjI+ib5tGjO3XinZxUbSPj6RdXBRxURraqvo00JWqhuyjx9m0P49N6Xls2p/H5v15/JBXdHJ/dEQI7eOj6N+mMdf3a02ip5edGBdJEw1t5Sca6EpV4UR4b07PY2Ml4d0hPoqBiU3okxBLr9axdGkeTePIUA1tVes00N1u3jynKwgoFcN70367nS+8e7aKIToi1MGKlTpFA93t4uOdrqDOyqlk2ORAJeHdu3UsvRM0vFXdp4HudjNn2seJE52swnEnwnvz/jw2plce3kka3irAaaC7XT0M9PyiEtbuPcRmz5DJpvRzh3ev1rH0bB1DjIa3cgENdOUK3+cUsjg1g+S0DFbtyj15FaWGt6pPNNBVQCorN6zfd5jk1AwWp2awPeMoAJ2bNeSeSztwWZd4erWO1fBW9YoGugoYhcWlfLUjm8VbM1i6LZPso8UEBwmDEpvw5LVtGdW9Ge3iopwuUynHaKCrOu2HvGMkp2aSnJrBN9/lUFxaTnRECJd3bcbI7s0Y3qUZsZHaC1cKNNDdb+FCpyu4IMYYthzIZ7FnKGXz/nwA2jaJ5LbB7RjVoxkDE5sQqku8KnUWDXS3i4x0uoIqFZWUsWJXDou3ZrAkLZMf8ooQgYvaNuax0d0Y1b0ZnZo11CsvlaqCBrrbvfqqfbz/fmfrOEP20eMsSbNDKV/tyKawuIzIsGAu7RzPw1d04fJuzYhvGO50mUoFFA10t5s71z46HOjGGHZkHrVDKVszWLfvMMZAi5gIbrioNSO7N2dIhzi9QbBSNaCBrvympKyc1btzWZSaQXJqJt/nFgLQu3UsvxzZmVHdm9OzVYwOpSjlIxroyucOFRQz/ctdvLdqL0eKSgkLCeKSjnH8fFgHRnZrTovYCKdLVMqVNNCVz+QdK+GNr3bx5jd7KCgu5ereLRnTtxWXdo7Xu8QrVQv0t0zV2NHjpbz19W5e+2oX+UWlXN27BQ+N6kKX5tFOl6ZUvaKB7nbLlvntSxcWl/LOir1M/+93HCosYVT35vzPFZ3p2SrWb20qpc5NA11dsKKSMt5b9T3Tlu0k+2gxw7o05eErutC3TSOnS1OqXtNAd7u//90+Pvpojb/U8dIy5q7exytLd5KRf5yhHeOYfnsXBrRrUuOvrZSqOa8CXURGAy8CwcDrxpg/n7G/LfA20MhzzOPGmMC65tytPv3UPtYg0EvKypm/Jp2Xl+xk/+FjDExszAs39WdIxzgfFamU8oUqA11EgoGpwBVAOrBaRBYYY7ZWOOw3wFxjzDQR6QEsBBL9UK+qRaVl5Xy8/gAvJu/g+9xC+rZpxJ9u6M2lneN17rhSdZA3PfRBwE5jzC4AEZkNjAUqBroBYjzPY4EDvixS1a7ycsMnG22Q78oqoGerGN64M4kR3ZppkCtVh3kT6K2BfRU+TgcGn3HMb4H/iMiDQBQwqrIvJCKTgEkAbdu2vdBalZ+Vlxu+2HKQ5xdvZ3vGUbo2j+Yft13Ej3u20CBXKgD46qTozcBMY8yzIjIEmCUivYwx5RUPMsbMAGYAJCUlGR+1rc6nQYMqDzHGkJyayXOLtrP1h3w6NI3ipZv7c23vlgQFaZArFSi8CfT9QJsKHyd4XqvobmA0gDFmhYhEAPFApi+KVDXw2Wfn3GWM4csd2Ty3aDsb9h2mbZNInv1pX8b2a0WIrjeuVMDxJtBXA51FpD02yCcAt5xxzPfASGCmiHQHIoAsXxaqfGv5d9k8v2g7q/cconWjBvz5ht6MG5CgN45QKoBVGejGmFIRmQJ8gZ2S+KYxZouIPA2kGGMWAI8Ar4nI/2BPkE40xuiQSl3wzDP28cknAUjZk8uz/9nOil05NI8J55mf9OKmpDaEhWiQKxXoxKncTUpKMikpKY60Xa8MHw7Ahln/4tlF2/lyexbxDcOYPLwTtw5uq+uPKxVgRGSNMSapsn16pajLbYlsxvMJl7B46jc0jgzlV1d14/Yh7XT1Q6VcSH+rXezfG3/gwd6307CsmEev7MLES9rTMFz/yZVyK/3tdqklaRn8cvY6ko7s57VtHxH79y+cLkkp5Wd6JsyFlu/M5r5319KjVQxvZC0jtlFDp0tSStUC7aG7zJq9udzzTgrt46J4+65BRE+Z7XRJSqlaoj10F9m8P4+Jb62meUwEs+4ZROOoMKdLUkrVIg10l9iRcYTb31hFTEQo790zmGbRnhsx/+pXdlNKuZ4OubjA3pwCbn19FaHBQbx3z2BaNaqwfsuKFc4VppSqVdpDD3AHDh/jltdWUVJWznv3DCYxPsrpkpRSDtFAD2CZR4q49fVV5B8rYdbdg+ncPNrpkpRSDtIhlwB1qKCY21//loN5Rbx7zyB6tY51uiSllMM00APQkaIS7nzrW3bnFPDWxIHnv0lzQkLtFaaUcpQGeoA5VlzG3TNT2Hogn+m3D+CSTvHn/4R3362dwpRSjtMx9AByvLSMSbNSSNmbywsT+jGye3OnS1JK1SHaQw8QJWXlTHl/HV/tyOav4/twbZ9W3n3iQw/Zxxde8F9xSqk6QQM9AJSVGx79cAOLtmbwuzE9uTGpTdWfdML69f4rTClVp+iQSx1njOHXH23i4/UHeGx0N+4cmuh0SUqpOkoDvQ4zxvDMp6nMXr2PB0d0YvLwjk6XpJSqwzTQ67DnF23nzW92c9cliTx8RReny1FK1XE6hl5HTVv2HS8t2cmEgW146toeiEj1vlAX/Y9AqfpCA70OemfFHv7yeRpj+rbiD9f3rn6YA8yY4bO6lFJ1mw651DEfpuzjqY+3cEWP5jx7Y1+Cg2oQ5kqpekUDvQ7598YfeGz+Ri7tHM8rt/QnNNgH/zyTJtlNKeV6OuRSR5y4qfOAdo2ZfvsAwkOCffOFt2/3zddRStV52kOvA067qfPEgUSG6f+zSqkLp4HusDNv6hwTEep0SUqpAKWB7iC9qbNSype8CnQRGS0i20Rkp4g8fo5jbhSRrSKyRUTe922Z7rMj4wh3vPktMRGhvFvxps6+1q+f3ZRSrlflYK2IBANTgSuAdGC1iCwwxmytcExn4FfAJcaYQyLSzF8Fu8GJmzoHBwnv3TOY1hVv6uxrusqiUvWGNz30QcBOY8wuY0wxMBsYe8Yx9wJTjTGHAIwxmb4t0z30ps5KKX/xJtBbA/sqfJzuea2iLkAXEflGRFaKyOjKvpCITBKRFBFJycrKql7FASzryHFuq3BT5y61cVPn226zm1LK9Xw1Py4E6AwMBxKAL0WktzHmcMWDjDEzgBkASUlJxkdtB4TDhcXc/sYqfqjtmzqnp9dOO0opx3nTQ98PVLyjQoLntYrSgQXGmBJjzG5gOzbgFZ6bOr/5LbuyC3j9zqTz39RZKaWqyZtAXw10FpH2IhIGTAAWnHHMv7C9c0QkHjsEs8uHdQa0X3ywji0H8nn1louqvqmzUkpVU5WBbowpBaYAXwCpwFxjzBYReVpExngO+wLIEZGtwFLg/xljcvxVdCDJzC9i6bYsHri8E6N66E2dlVL+49UYujFmIbDwjNeeqvDcAA97NlXB0m12ws9VvVs4U8CQIc60q5SqdbpoiJ8lp2bSKjaCrrUxo6Uyf/qTM+0qpWqdXvrvR8dLy/h6ZzYjujer2U0qlFLKCxrofrRqVy6FxWWM7Obg2Pm4cXZTSrmeDrn40ZK0TCJCgxjSMc65InL03LRS9YUGup8YY0hOy+CSjvFEhProZhUXXgQ0KIKyYPtch32UcjUNdD/5Luso+3KPcd+wjrXf+LHDsHEOrJkJg1Pta39sBY0ToUkHz2N7aNzefhzbBoL1R0GpQKe/xX6SnGqnK17etZYWnjQG0lfbEN/8Tyg9Bq36w/YEu/+OsXBoN+TshJ2LobTo1OcGhdhQb+IJ+MbtTwV+40QIi6yd96CUqhENdD9JTsuke8sYWvlzaVzw9Mbn2iDP3AJhDaHvBBgwEVr1g4xn7HGjnzz1OeXlcPQg5O6G3F026E88378GivJObyO65amQP9mz9zxG6jIGStUVGuh+kFdYwpq9h5jsr+EWYyA9Bda8dao33rIfXPci9BoH4RXmvD/55NmfHxQEMa3slnjJ2fsLcyuE/G7P812wM9n+R1BRRKPTh28qBn7DFrYtpVSt0ED3g//uyKKs3HB5Nx8PtxTlneqNZ2z29MZv8vTG+/uuncgmdms94Ox9xYVwaM/pPftDu+HAWtj6MZiyU8eGx0K3a+x/Mh2GQbDeL1Upf9JA94MlqRk0iQqjX5tGNf9ixthhkBO98ZJCaNkXrn0Beo8/vTdemauuso+ffVbzWsCOpzfvYbczlZVA3j4b9rm74cA6SP0UNrwPkXHQ4ye25jYXa89dKT/QQPexsnLDsu1ZjOjWjOCgGkwTPNkbfxsyNkFoFPT+qe2Nt77I+69z7Fj1a7hQwaGeYZcOp1679nnYsQg2z4f170PKGxDTGnpeb8O9ZT+dTqmUj2ig+9i67w9xuLCEEdUZbjEG9q+FNW+e6o236A3XPGfDPCLG9wX7W0g4dL/WbsePwrbPYPM8WDUdVrwCTTraYO81Dpp2dbpapQKaBrqPLUnLJCRIuLRzU+8/qSgfNnnGxg+e6I2P94yNX+SeHmx4Q+jzU7sV5kLqJzbc//tX+O9foHlv6D3Ohnujtk5Xq1TA0UD3sSVpmQxMbEJsgypOABpjTySumQmb5kNJQeD3xi9EZBMYcKfdjhyELR/ZYZnFv7VbwiD7n1rP66FhLc3lVyrAaaD70P7Dx0g7eIRfX9393AcV5cOmDz298Y0QGml7pAPusmPjvu6NX3utb7+eP0S3gIsn2+3QHhvsm+bDZ/8Lnz8O7S+DXuPtsE2Dxk5Xq1SdJfbeFLUvKSnJpKSkONK2v8xauZcn/7WZ5EeG0bFpw9N37l9rZ6qc6I0372WHVPrcCBG1dMPoQJOZ6gn3eXZqZHAYdBpl/wPsehWERTldoVK1TkTWGGOSKtunPXQfWpKaQWJcJB3izwiauXfYOdqhkdDrBk9vfIB7xsb9pVl3GPEbuPzXdnhq03zY8k/YttB+L7teZXvunUbak69K1XMa6D5yrLiM5d/lcMvgtqffzCIzzYb5wHth5JO13xsfPtw+LltWu+36koj9D7D1ALjy9/D9cttr3/qx7cFHxEL3Mbbn3v4yCHJodUulHKaB7iPLv8vmeGn52Tez2DgbJBiGPaZDK74QFASJP7Lb1X+DXctsuG/5CNbNgqhm9kRq3wkXNl9fKRfQQPeR5LRMosKCGdS+wmJV5WX24qBOo6DhBUxjVN4JDoXOV9it5Bhs/8JOg1wzE76dDm2HwJApdmhGe+2qHtBA9wFjDEvTMrm0c1PCQipc0r7nK8jfD1c+41xx9UVoA+j5E7sV5cG692DlNJhzq71y9eL7od8teiJVuZouqOEDqT8c4Ye8IkZ0P2O+9IY5EB4DXa92prD6KiIWhtwPv1gH49+yUx0XPgrP94TkZ+BIhtMVKuUX2kP3gSVpNiCGd60wrFJcYE/a9R5ne49OufFG59p2WnCInVXU83r4fqVdauCrZ2H5S9D7RhjyQOWLjCkVoDTQfSA5LZO+CbE0i4449WLqp3a+eZ8JzhUGcP/9zrZfF4hAuyF2y/kOVr5qh2TWvwsdR8LQKdDhcp1GqgJeQAb68JnDz3rtxp43cv/A+yksKeTq984e4pjYbyIT+00kuzCb8XPHn7V/ctJkbup1E/vy9nH7R7eftf+RIY9wXdfr2Ja9jZ9/+vOTr5eVNmDf9w8wdoD9Vq4/uJ6HPn+Ivx7cSZuQMG5Z+gRGhD+O/CND2wxl+b7lPJH8xFlf/4XRL9CvRT8W71rM77/8/Vn7p187na7xXflk2yc8u+LZs/bPun4WbWLbMGfzHKalTDv5evhxuz75e7d/RHxkPDPXz2Tm+plnff7CWxcSGRrJq6tfZe6WuWftXzZxGQB/X/53Pt3+6Wn7GoQ24LNb7fK8z/z3GZJ3J5+2Py4yjvk3zgfgV4t/xYr0FaftT4hJ4N0b3gXgoc8fYv3B9aft7xLXhRnXzQBg0ieT2J6z/bT9/Vr044XRLwBw2z9vIz0//bT9QxKG8KdRfwJgXPL/klOYQ0zLTow5ks31e/5L3HfJ9kKvIQ9w3ea3OVJ6/LTPv7bLtTw69FGgbv3snfCby37DqA6jTv7sncmpn70T5t04T3/2gHFzx5FTmHPae/K1gAz0uuRYQQdA6N321OmIuNISBhQd4b3Y5hiHe31/eX6TfXJ2TtRr+cEhvNuoBXNim/FggxaMyd0H/5rMzOBQ5kfH80l0PEf0xtkqwHh16b+IjAZeBIKB140xfz7HceOAecBAY8x5r+t3y6X/D7y3ltV7cln5q5EEnVj//JuXYNGTMGUNxHdytkA3XFhUG4yB75Jh+Suwa6m9ErX/bXZ9mYrruyvlsPNd+l/lLBcRCQamAlcBPYCbReSsM0kiEg38ElhVs3IDR0lZOV96bmZxMsyNgQ0fQMJA58NceU/EXi9wx7/gvm/s3ZVS3oKXLoI5t8H39ebHWgUwb6YtDgJ2GmN2GWOKgdnA2EqOewb4C1Dkw/rqtNV7cjlyvPT0e4ce3ASZW6HPTc4VpmqmRS+4fho8tAl+9D+w+yt480p4fRRs+Ze9YEypOsibQG8N7KvwcbrntZNE5CKgjTHm3+f7QiIySURSRCQlKyvrgouta5akZhIWHMSPOsWfenHjHAgKteuKqMAW0xJG/R88vBWu+hsUZMGHd8JL/WHlP+wdmJSqQ2p8YZGIBAHPAY9UdawxZoYxJskYk9S0aeBfCr9kWyYXd4wjKtxz8qys1F7q3+XH9gYOdcHEiXZT1RcWBYMnwYNr4cZZdv32zx+D53vAov+D/ANOV6gU4N0sl/1AmwofJ3heOyEa6AUs86wy2AJYICJjqjoxGsh2ZxewK6uAOy5ud+rFXUuhINMuDFVXaJj7TlAw9Bhjt32rYcXL9iKlFa/YZXyHTrF3nVLKId4E+mqgs4i0xwb5BOCWEzuNMXnAyTEHEVkGPOrmMAd7qzmAERVXV9ww215m3vlKh6qqRHa2fYyPP/9x6sK0GQht3rF3WFo5DdbOsitrth9mx907DNcLlVStq3LIxRhTCkwBvgBSgbnGmC0i8rSIjPF3gXXV0rRMOjdrSNu4SPtCUT6kfQo9b6hbN1sYP95uyj8aJ8JVf4GHt8Co30H2dpj1E3jransyVala5NWVE8aYhcDCM1576hzHDq95WXXb0eOlrNqdw88uaX/qxdQFUFoEfW92rjDlnAaN4UcP2Xnra9+xa8a8fS0kXgqXPwHthjpdoaoHdLXFavh6RxYlZYYRFacrbpgNTTpCQqXz/VV9ERIOg+61Kz2O/jNkbYO3roJ3fgL7vnW6OuVyGujVkJyaSUxECAPaee5Af/h7u/Z53wk6bqqs0Aa2t/7LDXDlH+z1CW9cAe+Og/Q1TlenXEoD/QKVlxuWbstkWNdmhAR7vn0bPQsK9anHS9WqyoVF2tkvD220Y+z718LrI+D9m+DA+qo/X6kLoKsPXaBN+/PIPlrMyBPDLcbY4Za2Q+0Jsrpm8mSnK1Bg57L/6CEYeDesmg7LX4YZw6DbtTD8cZ3uqHxCA/0CJadlEiQwrIvnwqgDayFnBwx90NnCzuUmXYKgTgmPhssetePsq6bbxcDSPoUeY2HY43rDDVUjOuRygZakZXBR28Y0jgqzL2yYDcHh9l6WddG+fXZTdUtELAz7XzsUM+wx2LkEpg2FD++yJ1KVqgYN9AuQkV/E5v35p+4dWloMm+ZBt6vtL2hddPvtdlN1U4NGdlrjQxvh0odh+xcwdTDMvxeydzpdnQowGugXYOnJq0M9gb5zMRzL1bnnquYim8DIp2ywX/ILOwwzdSB8NBlydzldnQoQGugXIDktk9aNGtC1ebR9YcMHENUUOo5wtjDlHlHxcMXTdrrjxffDln/Cy0nw8RQ4tNfp6lQdp4HupaKSMr7Zmc2Ibs0QETh2CLZ/bhdlCg51ujzlNg2bwY//YIN90L12auzLF8Env4TDek5EVU4D3UurdudSWFx2arhly0dQVly3VlZU7hPdwq4V88v1MOAuWP++XY/9349A3v6qP1/VKzpt0UtLUjOICA1iSMc4+8KG2dC0G7Ts62xhVXmkymXqVSCIaQXX/B0u+aVdJ2bNTLvCY9JddnXH6BZOV6jqAO2he8EYw5JtmfyoUzwRocH2JNW+VYFxqf9119lNuUOjNnDdC/ZmG31uhG9fgxf7wudPwNFMp6tTDtNA98LOzKPsyz126t6hG+YAAr0D4FL/bdvsptylcTsY+wo8mGKXbF41zQZ78tNQlOd0dcohGuheSK44XdEYz40MLoPY1lV8Zh3w85/bTblTkw72htYPrIauV9vhmBf72ZtulBY7XZ2qZRroXliSlkmPljG0jG1gh1oO7dG556puie8E49+AScvsujCfPw6vJNkL38rLna5O1RIN9CocLixmzd5Dp2a3bPgAQiOhu45LqzqoVX+442O4bT6Ex8D8u+G1y2HXMqcrU7VAA70K/92eRVm5sZf7lxTB5o9smIc3dLo0pSonAp1Gwc+/hOunQ2EOvDMWZt1g12VXrqWBXoWlaZnERYXRN6GRvZDoeJ7OPVeBISjI/qxOSYErnoH9a+Afl8JH9+nFSS6l89DPo6zcsGx7FiO6NSM4SOzc8+iW9s7ugeI3v3G6AuW00Ai7PsxFt8NXz9llezf/EwZPgksfsfdDVa6gPfTzWPf9IQ4XljCyW3MoyIadi6D3TyEo2OnSvDdqlN2UatAYrnwGHlwDvcbZtdhf7AvfvGiHE1XA00A/j+S0TEKChEu7xMPm+VBeGnizW9avt5tSJzRqY6c63vc1JAyERU/BywPssgLlZU5Xp2pAA/08lqRmMjCxCTERoXZ2S4vegXdHmYcesptSZ2rRy86GuWOBXeXxX5PtGPuORfZ6CxVwNNDPIf1QIdsyjjCyezN7B5kD6wKvd66UNzoMg3uXwvg3oaQA3hsPb19nb2itAooG+jmcdjOLDbNBguxSuUq5UVCQHVd/YDWM/gtkbrXz1+f9DHJ3O12d8pIG+jkkp2WSGBdJh7hIuxZ1x5EQ3dzpspTyr5AwuPg++MV6uPRRSFsIrwyEzx6zEwNUnaaBXonC4lKWf5fDiG7NYe/XkJ+uc89V/RIRAyOfhF+sg363wBCtdaMAAA8BSURBVLcz7BoxX/4Nigudrk6dg1eBLiKjRWSbiOwUkccr2f+wiGwVkY0ikiwi7Xxfau1ZvjOH4tLyU8MtYdHQ7Rqny6qeP/7RbkpVR0xLGPMSTF4B7S+FJb+3N9hYMxPKSp2uTp2hykAXkWBgKnAV0AO4WUTOnOqxDkgyxvQB5gF/9XWhtSk5LZOosGAGJUTA1o+h51gIbeB0WdUzdKjdlKqJZt3g5g/grs/ttMdPfgnThkLav3VGTB3iTQ99ELDTGLPLGFMMzAbGVjzAGLPUGHPi77CVQIJvy6w9xhiWpmVyWZemhO34DIqPBvbsluXL7aaUL7QbAncvghtngSmD2bfAW1fBvm+drkzhXaC3Biou/JDuee1c7gY+q2yHiEwSkRQRScnKyvK+ylq09Yd8DuYX2ZtZbPgAYttA2wDu4T7xhN2U8hUR6DEG7l8J1zwHOd/BG1fAh3fB4e+drq5e8+lJURG5DUgC/lbZfmPMDGNMkjEmqWnTpr5s2meWpNrpiiNbG9i1FPrcZKd0KaVOFxwKA++2J06HPQbbPoOXk2Dx7+D4Eaerq5e8Sar9QJsKHyd4XjuNiIwCfg2MMcYc9015tW/Jtkz6tmlE3O6PwZTr7BalqhLeEC5/wt4Or8dY+Po5eOkiWPO2LiVQy7wJ9NVAZxFpLyJhwARgQcUDRKQ/MB0b5gF7p9rso8dZv+8wI7p6Zre0HgDxnZ0uS6nAEJsA416De5ZA40T45BcwfRjs/tLpyuqNKgPdGFMKTAG+AFKBucaYLSLytIiM8Rz2N6Ah8KGIrBeRBef4cnXasm1ZGAPXNMuGjM2BfTJUKackDIC7/2OXEig6bJcR+OAWO9au/Mqr9dCNMQuBhWe89lSF565Yn3VpWibNY8Lp+MOnEBRi76Ye6F54wekKVH0kYpcS6Ho1rJgKXz8PUwfDoEkw7P/pGux+omf7PIpLy/lyexYjuzRBNn0InX8MUXFOl1Vz/frZTSknhDaAyx6FB9fa81ErX7Xj69++phcm+YEGukfKnlyOHC9lXOPv4GgG9L3J6ZJ8Y/FiuynlpOjmMPYVe5/T5j1h4aP2wqQdi5yuzFU00D2WpGUSFhJE35zPISIWuox2uiTf+P3v7aZUXdCyD9z5CUx4H8pL7FK9s26AzFSnK3MFDXSPJWmZDE9sQMj2f9uxv5Bwp0tSyp1E7NpI96+CK/8A6Skw7RL49yNQkON0dQFNAx3YnV3AruwC7ojdAKXHoI/OPVfK70LCYOgUe2FS0s8g5S278Nfyl6G02OnqApIGOrZ3DpCU9x9o3B7aDHK4IqXqkag4uObvMHk5tBkI//kNTB0EqZ/owl8XSAMdWJKWwdD4Y0Skf2Pnnos4XZJS9U+zbvYep7fOt0Oec26zc9h/2OB0ZQHDq3nobnakqIRvd+fyWvsUOGqgz41Ol+Rb06c7XYFSF6bzKOgwHNa8BUv/aK827X8rjHgSols4XV2dVu976F/vyKakrJzBRxZB2yHQpL3TJflW1652UyqQBIfAoHvt+PqQB2DDHDt//cu/Qckxp6urs+p9oCenZTIkYh8N8nbalRXd5pNP7KZUIGrQCH78B3hgFXS83N4x6ZWBsGmejq9Xol4Henm5Ydm2TO5r9C0Eh0PPnzhdku89+6zdlApkcR1hwnt2DnuDRjD/bnjjSjvlUZ1UrwN94/48Dh8t5OLCpdD1Kl1fQqm6rv1lMOm/MOZlOLQHXh8J8++F/ANOV1Yn1OtAX5KawfDgDYQXH9J1z5UKFEHBcNEd8Iu1cOkj9r6/LyfBl3+HkiKnq3NU/Q70bZn8LPpbiIyDTq5YMFKp+iM8GkY+VWF8/Rl4dTCkLay34+v1NtAP5hXx/f4DDDq+CnqNt7fTUkoFnibt7fj67R/Zc2Gzb4Z3b4CsbU5XVuvqbaAv3ZbJNcGrCDHF7h5umTXLbkq5XccRMPkbGP1nSF9jV3P8/AkoynO6slpTbwN9SVomE8KWY+K7Qqv+TpfjP23a2E2p+iA4FC6ebMfX+916av31te9AebnT1fldvQz0opIy9u7YTF+TivS9yd2X+s+ZYzel6pOoeBjzEkxaaqc8LngQXrscvl/ldGV+VS8DfeWuHK4q/xKDQG+XXep/pmnT7KZUfdSqP/zsC7jhdXvjmjevhH9OgvwfnK7ML+ploC9NzeCGkK8pT7wUGulwhFKuJgJ9fgpTUuw0xy0fwcsD4KvnoPS409X5VL0LdGMMGVu/op1kEOzmk6FKqdOFNzw1zbHDcEj+nb1x9bbPXDPNsd4F+o7Mo/yocDGlQRHQY4zT5SilaluTDnDz+3DbP+1J1A8mwLvjIGu705XVWL0L9GVb07kueAUlXa62FyYopeqnTiPtTTV+/CdIXw3ThsAXvw7oaY71LtCPbPiUWCmkQdKtTpdSO+bNs5tS6mzBoTDkfnhwrb25zYqpdnx97ayAnOZYrwL9cGExfXI/52hoPLQf7nQ5tSM+3m5KqXNr2BTGvmKnOTZuDwum2IW/9q12urILUq8CfcWm7QyXdRR0vd4uoF8fzJxpN6VU1Vr1h7v/A9fPsCs4vjEKProPjhx0ujKv1KtAL1g7h1Apo+kldzpdSu3RQFfqwohA35vgwRT40f/A5vl2GObr5+v8NEevAl1ERovINhHZKSKPV7I/XETmePavEpFEXxdaU6Vl5XTNWMj+8E4EteztdDlKqbouPBpG/RbuXwmJl8Li38KrF8O2z+vsNMcqA11EgoGpwFVAD+BmEelxxmF3A4eMMZ2A54G/+LrQmkrdvIbe7CS/yw1Ol6KUCiRxHeGW2XDrfJBg+OAmeO+nkL3D6crO4k0PfRCw0xizyxhTDMwGxp5xzFjgbc/zecBIkbq1QEr+qvcoM0LCsDucLkUpFYg6j7LTHK/8A+xbZXvrXz/vdFWn8SbQWwP7Knyc7nmt0mOMMaVAHhB35hcSkUkikiIiKVlZWdWruJpCYpqzJu5aouP1Un+lVDWFhMHQKfDgGrvsdkyC0xWdplanehhjZgAzAJKSkmp1EGrwhF/VZnN1x8KFTleglPs0bAZjpzpdxVm86aHvByp2axM8r1V6jIiEALFAji8KVDUUGWk3pZTreRPoq4HOItJeRMKACcCCM45ZAJyYCzgeWGJMHT0NXN+8+qrdlFKuV2Wge8bEpwBfAKnAXGPMFhF5WkROrG71BhAnIjuBh4GzpjYqh8ydazellOt5NYZujFkILDzjtacqPC8Cfurb0pRSSl2IenWlqFJKuZkGulJKuYQGulJKuUQ9WXKwHlu2zOkKlFK1RHvoSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEhroSinlEuLUsuUikgXsreVm44Fsl7bn5vdW2+25+b3Vdntufm9OtAfQzhjTtLIdjgW6E0QkxRiT5Mb23Pzears9N7+32m7Pze/NifaqokMuSinlEhroSinlEvUt0Ge4uD03v7fabs/N762223Pze3OivfOqV2PoSinlZvWth66UUq6lga6UUi7hykAXkdEisk1EdorI45XsDxeROZ79q0Qk0Y9tXSYia0WkVETGV7edC2jvYRHZKiIbRSRZRNr5ub37RGSTiKwXka9FpIc/26tw3DgRMSJS7SljXry3iSKS5Xlv60Xknuq25U17nmNu9Pz7bRGR9/3Znog8X+G9bReRw35sq62ILBWRdZ6fzaur25aX7bXz/PxvFJFlIpJQg7beFJFMEdl8jv0iIi95atkoIhdVt60aM8a4agOCge+ADkAYsAHoccYx9wP/8DyfAMzxY1uJQB/gHWB8Lby3y4FIz/PJ1X1vF9BeTIXnY4DP/dme57ho4EtgJZDkx/c2EXilFn8uOwPrgMaej5v5+3tZ4fgHgTf9+N5mAJM9z3sAe/z8vfwQuNPzfAQwqwbtXQZcBGw+x/6rgc8AAS4GVvniZ6Y6mxt76IOAncaYXcaYYmA2MPaMY8YCb3uezwNGioj4oy1jzB5jzEagvBpfvzrtLTXGFHo+XAlUu2fiZXv5FT6MAmpylt2bfzuAZ4C/AEW10JaveNPevcBUY8whAGNMpp/bq+hm4AM/tmWAGM/zWOBANdvytr0ewBLP86WV7PeaMeZLIPc8h4wF3jHWSqCRiLSsbns14cZAbw3sq/Bxuue1So8xxpQCeUCcn9rypQtt725sz8Gv7YnIAyLyHfBX4Bf+bM/z52wbY8y/a9COV215jPP8GT1PRNr4ub0uQBcR+UZEVorIaD+3B9jhCaA9pwLQH239FrhNRNKBhdi/CKrLm/Y2ADd4nl8PRItIdX7HfVVPrXBjoCtARG4DkoC/+bstY8xUY0xH4DHgN/5qR0SCgOeAR/zVxhk+ARKNMX2ARZz6q85fQrDDLsOxPebXRKSRn9sEO+w4zxhT5sc2bgZmGmMSsEMUszz/nv7yKDBMRNYBw4D9gD/fX53gxkDfD1TsSSV4Xqv0GBEJwf4JmOOntnzJq/ZEZBTwa2CMMea4v9urYDbwEz+2Fw30ApaJyB7seOWCap4YrfK9GWNyKnz/XgcGVKMdr9vD9uwWGGNKjDG7ge3YgPdXeydMoPrDLd62dTcwF8AYswKIwC5s5Zf2jDEHjDE3GGP6Y38XMMZU+6RvTeupNU4N3vtrw/ZydmH/hDxxwqTnGcc8wOknRef6q60Kx86k5idFvXlv/bEnjDrX0veyc4Xn1wEp/mzvjOOXUf2Tot68t5YVnl8PrPTz93I08LbneTz2z/g4f34vgW7AHjwXGfrxvX0GTPQ8744dQ69Wm162Fw8EeZ7/AXi6hr8LiZz7pOg1nH5S9NuatFWjOp1q2K9vyv5Jt90TbL/2vPY0tscKtnfwIbAT+Bbo4Me2BmJ7XgXYvwK2+Pm9LQYygPWebYGf23sR2OJpa2lloeHL9s44dhnVDHQv39ufPO9tg+e9dfPz91KwQ0pbgU3ABH9/L7Fj23+uSTtevrcewDee7+V64Eo/tzce2OE55nUgvAZtfQD8AJR4fpfvBu4D7qvw7zbVU8ummvxM1nTTS/+VUsol3DiGrpRS9ZIGulJKuYQGulJKuYQGulJKuYQGulJKuYQGulJKuYQGulJKucT/B+MR5Uj6vXl0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xticks(np.arange(0,1.1,0.1))\n",
    "plt.plot(x_normalized, y_normalized);\n",
    "plt.plot(x_difference, y_difference);\n",
    "plt.hlines(Tmx[0], plt.xlim()[0], plt.xlim()[1], colors='g', linestyles='dashed');\n",
    "plt.vlines(x_difference_maxima, plt.ylim()[0], plt.ylim()[1], colors='r', linestyles='dashed');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertical, red dashed line represents the x value of the knee point. The horizontal greeb dashed line represents the threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2222222222222222"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "knee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2222222222222222"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# normalized x value where the knee was determined\n",
    "norm_knee"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "kneed",
   "language": "python",
   "name": "kneed"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
   }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}